A start on v.16

.gitignore

@@ -0,0 +1,5 @@
+.Rproj.user
+.Rhistory
+.RData
+.Rbuildignore
+demography.Rproj

ChangeLog

@@ -1,3 +1,8 @@
+v1.16
+- Fixed bug in pop.sim when no migration data was used
+- Added hmd.pop() to read population data from www.mortality.org.
+- Fixed a bug in forecast.fdm() when the time series frequency is greater than 1.
+
 v1.15
 - smooth.demogdata will no longer return NAs for fertility data. Instead, the fertility rate for the nearest age with positive rate is used.
 - Fixed occasional bug in computing life expectancy prediction intervals from coherent fdm model.

DESCRIPTION

@@ -1,6 +1,6 @@
 Package: demography
-Version: 1.15
-Date: 2013-08-20
+Version: 1.16
+Date: 2013-08-28
 Title: Forecasting mortality, fertility, migration and population data
 Description: Functions for demographic analysis including lifetable
         calculations; Lee-Carter modelling; functional data analysis of

NAMESPACE

@@ -1,7 +1,7 @@
 export(bms, cm.spline, cm.splinefun, coherentfdm, combine.demogdata,
 compare.demogdata, demogdata, extract.ages, extract.years, fdm,
 fitted.fdm, fitted.lca, forecast.fdm, forecast.lca, forecast.fdmpr,
-isfe, isfe.demogdata, lca, hmd.mx,hmd.e0,
+isfe, isfe.demogdata, lca, hmd.mx,hmd.e0,hmd.pop,
 life.expectancy, lifetable, mean.demogdata,
 median.demogdata, netmigration, plot.demogdata, plot.errorfdm,
 plot.lifetable, read.demogdata,flife.expectancy,

R/demogdata.R

@@ -22,7 +22,7 @@ demogdata <- function(data, pop, ages, years, type, label, name, lambda)
     if(length(years) != n)
         stop("Number of years doesn't match data")
 
-    types <- c("mortality","fertility","migration")
+    types <- c("mortality","fertility","migration","population")
     idx <- pmatch(type,types)
     if(is.na(idx))
         warning("Unknown type")
@@ -39,10 +39,19 @@ demogdata <- function(data, pop, ages, years, type, label, name, lambda)
             lambda <- 1
     }
 
-    obj <- list(year=years, age=ages, rate=list(as.matrix(data)), pop=list(as.matrix(pop)), type=type,
-        label=label, lambda=lambda)
-    dimnames(obj$rate[[1]]) <- dimnames(obj$pop[[1]]) <- list(ages,years)
-    names(obj$rate) <- names(obj$pop) <- name
+    if(type=="population")
+    {
+        obj <- list(year=years, age=ages, pop=list(as.matrix(pop)), type=type,
+            label=label, lambda=lambda)
+        dimnames(obj$pop[[1]]) <- list(ages,years)
+    }
+    else
+    {
+        obj <- list(year=years, age=ages, rate=list(as.matrix(data)), pop=list(as.matrix(pop)), type=type,
+            label=label, lambda=lambda)
+        dimnames(obj$rate[[1]]) <- dimnames(obj$pop[[1]]) <- list(ages,years)
+        names(obj$rate) <- names(obj$pop) <- name
+    }
     return(structure(obj,class="demogdata"))
 }
 
@@ -138,6 +147,8 @@ plot.demogdata <- function(x, series=ifelse(!is.null(x$rate),names(x$rate)[1],na
     series <- tolower(series)
     ages <- ages[ages <= max.age]
     data <- extract.ages(extract.years(x,years),ages,FALSE)
+    if(x$type == "population")
+        datatype <- "pop"
 
     # Extract data matrix
     if(!is.element(datatype,names(data)))
@@ -226,6 +237,9 @@ lines.demogdata <- function(x, series=ifelse(!is.null(x$rate),names(x$rate)[1],n
     ages <- ages[ages <= max.age]
     data <- extract.ages(extract.years(x,years),ages,FALSE)
 
+    if(x$type == "population")
+        datatype <- "pop"
+
     # Extract data matrix
     if(!is.element(datatype,names(data)))
         stop(paste("Data type",datatype,"not found"))

R/fdm.R

@@ -178,7 +178,7 @@ forecast.fdm <- function(object,h=50,level=80, jumpchoice=c("fit","actual"),
     output <- list(
         label=object$label,
         age=object$age,
-        year=max(object$year)+(1:h),
+        year=max(object$year)+(1:h)/tsp(object$year)[3],
         rate=list(forecast=fcast$mean$y,
                   lower=fcast$lower$y,
                   upper=fcast$upper$y),

R/hmd.R

@@ -59,3 +59,38 @@ hmd.e0 <- function(country, username, password)
     return(lt)
 }
 
+
+hmd.pop <- function(country, username, password, label=country)
+{
+    path <- paste("http://www.mortality.org/hmd/", country, "/STATS/", "Population.txt", sep = "")
+    userpwd <- paste(username, ":", password, sep = "")
+    txt <- RCurl::getURL(path, userpwd = userpwd)
+    con <- textConnection(txt)
+    pop <- try(read.table(con, skip = 2, header = TRUE, na.strings = "."),TRUE)
+    close(con)
+    if(class(pop)=="try-error")
+        stop("Population file not found at www.mortality.org")
+
+    obj <- list(type="population",label=label,lambda=0)
+
+    obj$year = sort(unique(pop[, 1]))
+    n <- length(obj$year)
+    m <- length(unique(pop[, 2]))
+    obj$age <- pop[1:m, 2]
+    mnames <- names(pop)[-c(1, 2)]
+    n.pop <- length(mnames)
+    obj$pop <- list()
+    for (i in 1:n.pop) 
+    {
+        obj$pop[[i]] <- matrix(pop[, i + 2], nrow = m, ncol = n)
+        obj$pop[[i]][obj$pop[[i]] < 0] <- NA
+        dimnames(obj$pop[[i]]) <- list(obj$age, obj$year)
+    }
+    names(obj$pop) <- tolower(mnames)
+
+    obj$age <- as.numeric(as.character(obj$age))
+    if (is.na(obj$age[m])) 
+        obj$age[m] <- 2 * obj$age[m - 1] - obj$age[m - 2]
+    return(structure(obj, class = "demogdata"))
+}
+

R/lca.R

@@ -118,14 +118,14 @@ lca <-  function(data,series=names(data$rate)[1],years=data$year, ages=data$age,
     else if(adjust=="e0")
     {
         e0 <- apply(mx,1,get.e0,agegroup=agegroup,sex=series,startage=startage)
-        FUN <- function(p,e0i,ax,bx,agegroup,series,startage){e0i - estimate.e0(p,ax,bx,agegroup,series,startage)}
+        FUN2 <- function(p,e0i,ax,bx,agegroup,series,startage){e0i - estimate.e0(p,ax,bx,agegroup,series,startage)}
         for (i in 1:m)
         {
             if(i==1)
                 guess <- kt[1]
             else
                 guess <- mean(c(ktadj[i-1],kt[i]))
-            ktadj[i] <- findroot(FUN, guess=guess, margin = 3*ktse[i],e0i=e0[i],ax=ax,bx=bx,agegroup=agegroup,series=series,startage=startage)
+            ktadj[i] <- findroot(FUN2, guess=guess, margin = 3*ktse[i],e0i=e0[i],ax=ax,bx=bx,agegroup=agegroup,series=series,startage=startage)
         }
     }
     else if(adjust=="none")