Nothing
#' @export
netmigration <- function(mort, fert, startyearpop=mort, mfratio = 1.05)
{
# Basic checks on inputs
if (!inherits(mort, "demogdata") | !inherits(fert, "demogdata"))
stop("Inputs not demogdata objects")
if (mort$type != "mortality")
stop("mort not mortality data")
if (fert$type != "fertility")
stop("fert not fertility data")
# Find years with both mortality and fertility data and startpop data
yrs <- mort$year[sort(match(fert$year, mort$year))]
yrs <- yrs[sort(match(startyearpop$year,yrs))]
if(max(startyearpop$year) > max(yrs))
startyearpop <- extract.years(startyearpop, c(yrs,max(yrs)+1))
else
{
startyearpop <- extract.years(startyearpop, yrs)
yrs <- sort(unique(pmin(yrs, max(yrs)-1)))
}
mort <- extract.years(mort, yrs)
fert <- extract.years(fert, yrs)
n <- length(yrs)
p <- length(mort$age)
# Splits births by male and female
B <- colSums(fert$pop$female * fert$rate$female/1000)
Bf <- B * 1/(1 + mfratio)
Bm <- B * mfratio/(1 + mfratio)
# Compute non-survival ratios from mortality rates
nsr.f <- 1 - lifetable(mort, "female", max.age = max(mort$age))$rx
nsr.m <- 1 - lifetable(mort, "male", max.age = max(mort$age))$rx
# Estimate deaths
oldest.f <- colSums(as.matrix(startyearpop$pop$female[(p - 1):p, ]))
oldest.m <- colSums(as.matrix(startyearpop$pop$male[(p - 1):p, ]))
if (n > 1)
{
Df <- nsr.f * rbind(Bf, startyearpop$pop$female[1:(p - 2), -n - 1], oldest.f[-n - 1])
Dm <- nsr.m * rbind(Bm, startyearpop$pop$male[1:(p - 2), -n - 1], oldest.m[-n - 1])
} else
{
Df <- nsr.f * c(Bf, startyearpop$pop$female[1:(p - 2), -n - 1], oldest.f[-n - 1])
Dm <- nsr.m * c(Bm, startyearpop$pop$male[1:(p - 2), -n - 1], oldest.m[-n - 1])
}
Dm[is.na(Dm) | Dm < 0] <- 0
Df[is.na(Df) | Df < 0] <- 0
# Compute net migration
Mf <- Mm <- matrix(NA, nrow = p, ncol = n)
for (j in 1:n)
{
current <- extract.years(startyearpop, years = yrs[j] + 1)
prev <- extract.years(startyearpop, years = yrs[j])
Mf[, j] <- current$pop$female - c(Bf[j], prev$pop$female[1:(p - 2), ], oldest.f[j]) + Df[, j]
Mm[, j] <- current$pop$male - c(Bm[j], prev$pop$male[1:(p - 2), ], oldest.m[j]) + Dm[, j]
}
# Store migration figures in a demogdata object with same population figures as startyearpop
mig <- extract.years(startyearpop, years = yrs)
mig$rate$female <- Mf
mig$rate$male <- Mm
mig$rate$total <- Mm + Mf
mig$pop$total <- mig$pop$male + mig$pop$female
mig$lambda <- 1
dimnames(mig$rate$male) <- dimnames(mig$rate$female) <- dimnames(mig$rate$total) <- dimnames(mig$pop$male)
mig$type <- "migration"
# Return result. Note: rows are actually (B,0), (0,1), etc.
return(mig)
}
# data must be a demogdata object containing population values for the last year of observation mort
# and mig are fmforecast2 objects and fert is a fmforecast object. If they are NULL, it is assumed
# all values are zero.
#' Population simulation
#'
#' Simulate future sample paths of a population using functional models for mortality, fertility and migration.
#'
#' @param mort Forecasts of class \code{fmforecast2} for mortality.
#' @param fert Forecasts of class \code{fmforecast} for female fertility.
#' @param mig Forecasts of class \code{fmforecast2} for net migration.
#' @param firstyearpop Population for first year of simulation.
#' @param N Number of sample paths to simulate.
#' @param mfratio Male-female ratio used in distributing births.
#' @param bootstrap If TRUE, simulation uses resampled errors rather than normally distributed errors.
#'
#' @return A list of two arrays containing male and female future simulated population values.
#' The arrays are of dimension (p,h,N) where p is the number of age groups, h is the forecast horizon
#' and N is the number of simulated sample paths.
#' @author Rob J Hyndman
#' @seealso \code{\link{simulate.fmforecast}}, \code{\link{simulate.fmforecast2}}.
#' @examples
#' \dontrun{
#' require(addb)
#' # Construct data objects
#' mort.sm <- smooth.demogdata(set.upperage(extract.years(australia,1950:2002),100))
#' fert.sm <- smooth.demogdata(extract.years(aus.fertility,1950:2002))
#' aus.mig <- netmigration(set.upperage(australia,100),aus.fertility,mfratio=1.0545)
#' # Fit models
#' mort.fit <- coherentfdm(mort.sm)
#' fert.fit <- fdm(fert.sm)
#' mig.fit <- coherentfdm(aus.mig)
#' # Produce forecasts
#' mort.fcast <- forecast(mort.fit)
#' fert.fcast <- forecast(fert.fit)
#' mig.fcast <- forecast(mig.fit)
#' # Simulate
#' aus.sim <- pop.sim(mort.fcast,fert.fcast,mig.fcast,australia)}
#'
#' @keywords models
#' @export
pop.sim <- function(mort, fert = NULL, mig = NULL, firstyearpop, N = 100, mfratio = 1.05, bootstrap = FALSE)
{
no.mortality <- FALSE # Not possible to proceed without mort object
no.fertility <- is.null(fert)
no.migration <- is.null(mig)
# Basic checks on inputs
if (!no.mortality)
{
if (!inherits(mort, "fmforecast2"))
stop("Inputs not fmforecast2 objects")
if (mort$female$type != "mortality" | mort$male$type != "mortality")
stop("mort not based on mortality data")
}
if (!no.fertility)
{
if (!inherits(fert, "fmforecast"))
stop("Inputs not fmforecast objects")
if (fert$type != "fertility")
stop("fert not based on fertility data")
}
if (!no.migration)
{
if (!inherits(mig, "fmforecast2"))
stop("Inputs not fmforecast2 objects")
if (mig$male$type != "migration" | mig$female$type != "migration")
stop("mig not based on migration data")
}
firstyr <- mort$male$year
if (!no.fertility)
firstyr <- intersect(firstyr, fert$year)
if (!no.migration)
firstyr <- intersect(firstyr, mig$male$year)
firstyr <- min(firstyr)
pop <- extract.years(firstyearpop, firstyr)$pop
# Check ages match. First make them all integers to prevent integer/numeric clashes
firstyearpop$age <- as.integer(firstyearpop$age)
mort$male$age <- as.integer(mort$male$age)
mig$male$age <- as.integer(mig$male$age)
if (!no.migration)
{
if (!identical(firstyearpop$age, mig$male$age))
stop("Please ensure that migration and population data have the same age dimension")
}
if (!no.mortality)
{
if (!identical(firstyearpop$age, mort$male$age))
stop("Please ensure that mortality and population data have the same age dimension")
}
p <- length(firstyearpop$age)
# Simulate all components
hm <- hf <- h <- Inf
if (!no.mortality)
{
mort.sim <- simulate(mort, nsim = N, bootstrap = bootstrap)
hm <- length(mort$male$year)
}
if (!no.fertility)
{
fert.sim <- simulate(fert, nsim = N, bootstrap = bootstrap)
hf <- length(fert$year)
}
if (!no.migration)
{
mig.sim <- simulate(mig, nsim = N, bootstrap = bootstrap)
h <- length(mig$male$year)
nm <- length(mig$male$model$year)
}
h <- min(hm, hf, h)
# Set up storage space
if (!no.fertility)
fage <- is.element(rownames(pop$female), fert$age)
pop.f <- pop.m <- array(0, c(p, h, N))
dimnames(pop.f) <- dimnames(pop.m) <- list(mort$female$age, 1:h, 1:N)
advance <- function(x0, x)
{
n <- length(x)
newx <- c(x0, x[1:(n - 2)], x[n - 1] + x[n])
}
# Simulate N future sample paths of population numbers
for (i in 1:N)
{
# Start with final observed populations
popf <- round(c(pop$female))
popm <- round(c(pop$male))
for (j in 1:h)
{
# Compute net migration
if (no.migration)
{
netf <- netm <- 0
Rf <- popf
Rm <- popm
} else
{
# Simulate net migration
netf <- round(mig.sim$female[, j, i] + mig$female$model$res$y[, sample(1:nm, 1)])
netm <- round(mig.sim$male[, j, i] + mig$male$model$res$y[, sample(1:nm, 1)])
# Add half migrants to current population
Rf <- pmax(popf + 0.5 * c(netf[2:(p - 1)], 0.5 * netf[p], 0.5 * netf[p]), 0)
Rm <- pmax(popm + 0.5 * c(netm[2:(p - 1)], 0.5 * netm[p], 0.5 * netm[p]), 0)
}
# Survivorship ratios
# firstyear pop used only to get structure. Data replaced.
Rt <- firstyearpop
Rt$type <- "mortality"
Rt$lambda <- 0
Rt$year <- 1
Rt$rate$female <- matrix(mort.sim$female[, j, i], ncol = 1)
Rt$rate$male <- matrix(mort.sim$male[, j, i], ncol = 1)
Rt$pop$female <- matrix(Rf, ncol = 1)
Rt$pop$male <- matrix(Rm, ncol = 1)
Rt$rate$total <- Rt$pop$total <- NULL
colnames(Rt$pop$male) <- colnames(Rt$pop$female) <- colnames(Rt$rate$male) <- colnames(Rt$rate$female) <- "1"
rownames(Rt$pop$male) <- rownames(Rt$pop$female) <- rownames(Rt$rate$male) <- rownames(Rt$rate$female) <- Rt$age
nsr.f <- 1 - lifetable(Rt, "female", max.age = max(Rt$age))$rx
nsr.m <- 1 - lifetable(Rt, "male", max.age = max(Rt$age))$rx
# Simulate deaths
cohDf <- pmax(nsr.f[c(2:p, p)] * Rf, 0)
cohDm <- pmax(nsr.m[c(2:p, p)] * Rm, 0)
Rf2 <- advance(0, Rf - cohDf) # Ignore deaths to births for now
Rm2 <- advance(0, Rm - cohDm)
Ef <- 0.5 * (Rf + Rf2)
Em <- 0.5 * (Rm + Rm2)
Df <- rpois(rep(1, p), Ef * mort.sim$female[, j, i])
Dm <- rpois(rep(1, p), Em * mort.sim$male[, j, i])
# Compute adjusted population ignoring births
cohDf[2:(p - 2)] <- 0.5 * (Df[2:(p - 2)] + Df[3:(p - 1)])
cohDm[2:(p - 2)] <- 0.5 * (Dm[2:(p - 2)] + Dm[3:(p - 1)])
cohDf[p - 1] <- 0.5 * Df[p - 1] + Df[p]
cohDm[p - 1] <- 0.5 * Dm[p - 1] + Dm[p]
cohDf[p] <- cohDm[p] <- 0
Rf2 <- advance(0, Rf - cohDf) # Fix problem with deaths to births later
Rm2 <- advance(0, Rm - cohDm)
# Simulate births from Poisson distribution
if (no.fertility)
B <- 0
else
{
lambda <- 0.5 * (Rf[fage] + Rf2[fage]) * fert.sim[, j, i]/1000
B <- sum(rpois(rep(1, length(fert$age)), lambda))
}
# Randomly split births into two sexes
Bm <- rbinom(1, B, mfratio/(1 + mfratio))
Bf <- B - Bm
# Infant mortality
RfB <- pmax(Bf + 0.5 * netf[1], 0)
RmB <- pmax(Bm + 0.5 * netm[1], 0)
cohDfB <- pmax(nsr.f[1] * RfB, 0)
cohDmB <- pmax(nsr.m[1] * RmB, 0)
Rf20 <- RfB - cohDfB
Rm20 <- RmB - cohDmB
Ef0 <- 0.5 * (Rf[1] + Rf20)
Em0 <- 0.5 * (Rm[1] + Rm20)
Df0 <- rpois(1, Ef0 * mort.sim$female[1, j, i])
Dm0 <- rpois(1, Em0 * mort.sim$male[1, j, i])
f0f <- cohDfB/(Ef0 * mort.sim$female[1, j, i])
f0m <- cohDmB/(Em0 * mort.sim$male[1, j, i])
cohDfB <- f0f * Df0
cohDmB <- f0m * Dm0
# Now we can fix infant deaths
cohDf[1] <- (1 - f0f) * Df0 + 0.5 * Df[1]
cohDm[1] <- (1 - f0m) * Dm0 + 0.5 * Dm[1]
Rf2 <- advance(RfB - cohDfB, Rf - cohDf)
Rm2 <- advance(RmB - cohDmB, Rm - cohDm)
# Add remaining migrants
popf <- round(pmax(Rf2 + 0.5 * netf, 0))
popm <- round(pmax(Rm2 + 0.5 * netm, 0))
# Store final population numbers
pop.f[, j, i] <- popf
pop.m[, j, i] <- popm
}
}
# Set up names of years
dimnames(pop.m)[[2]] <- dimnames(pop.m)[[2]] <- firstyr + (1:h)-1
# Return all sample paths
return(list(male = pop.m, female = pop.f))
}
deaths <- function(x)
{
if (!inherits(x, "demogdata") | x$type != "mortality")
stop("Not a mortality object")
npop <- length(x$rate)
out <- list()
for (i in 1:npop)
{
out[[i]] <- round(x$rate[[i]] * x$pop[[i]])
out[[i]][is.na(out[[i]])] <- 0
}
names(out) <- names(x$rate)
return(out)
}
births <- function(x)
{
if (!inherits(x, "demogdata") | x$type != "fertility")
stop("Not a fertility object")
out <- round(x$rate[[1]] * x$pop[[1]]/1000)
out[is.na(out)] <- 0
return(out)
}
Any scripts or data that you put into this service are public.
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.