Nearly finished v1.15

ChangeLog

@@ -1,5 +1,8 @@
 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.
+- Changed the way missing values are handled at the ends of the age range when smoothing.
+- Allowed missing values when using fdm().
 
 v1.14
 - minor changes to lifetable calculation.

DESCRIPTION

@@ -1,16 +1,16 @@
 Package: demography
-Version: 1.16
-Date: 2013-01-??
+Version: 1.15
+Date: 2013-08-20
 Title: Forecasting mortality, fertility, migration and population data
 Description: Functions for demographic analysis including lifetable
         calculations; Lee-Carter modelling; functional data analysis of
         mortality rates, fertility rates, net migration numbers; and
         stochastic population forecasting.
-Depends: R (>= 2.15.2), forecast (>= 3.09), rainbow, ftsa
-Imports: mgcv, cobs, strucchange, RCurl
+Depends: R (>= 2.15.2), forecast (>= 3.09), rainbow, ftsa, cobs
+Imports: mgcv, strucchange, RCurl
 LazyData: yes
 ByteCompile: TRUE
 Author: Rob J Hyndman with contributions from Heather Booth, Leonie Tickle and John Maindonald.
 Maintainer: Rob J Hyndman <[email protected]>
 License: GPL (>= 2)
-URL: http://robjhyndman.com/software/demography/
+URL: http://robjhyndman.com/software/demography/

NAMESPACE

@@ -49,3 +49,5 @@ import(forecast)
 import(ftsa)
 importFrom("stats", "simulate")
 importFrom("stats", "update")
+import(mgcv)
+import(cobs)

R/coherent.R

@@ -5,33 +5,33 @@ coherentfdm <- function(data, order1=6, order2=6, ...)
 {
 	# Check if missing data
 
-    # Use male and female if available
-    gps <- names(data$rate)
-    if(is.element("male",gps) & is.element("female",gps))
-    {
-        notneeded <- (1:length(gps))[-match(c("male","female"),gps)]
-        for(i in notneeded)
-            data$rate[[i]] <- data$pop[[i]] <- NULL
-    }
+  # Use male and female if available
+  gps <- names(data$rate)
+  if(is.element("male",gps) & is.element("female",gps))
+  {
+    notneeded <- (1:length(gps))[-match(c("male","female"),gps)]
+    for(i in notneeded)
+      data$rate[[i]] <- data$pop[[i]] <- NULL
+  }
+
+  J <- length(data$rate)
+  rate.ratio <- fdm.ratio <- list()
 
-    J <- length(data$rate)
-    rate.ratio <- fdm.ratio <- list()
-
-    # Construct ratio and product objects
+  # Construct ratio and product objects
 	is.mortality <- (data$type=="mortality")
-    y <- as.numeric(is.mortality)
-    pop.total <- 0
-    for (j in 1:J) 
-    {
-		if(is.mortality)
-			y <- y * data$rate[[j]]^(1/J)
-		else
-			y <- y + data$rate[[j]]/J
-        pop.total <- pop.total + data$pop[[j]]
-    }
-    rate.product <- demogdata(y, pop=pop.total, data$age, data$year, type=data$type, 
-        label=data$label, name="product")
-    for (j in 1:J)
+  y <- as.numeric(is.mortality)
+  pop.total <- 0
+  for (j in 1:J) 
+  {
+  	if(is.mortality)
+  		y <- y * data$rate[[j]]^(1/J)
+  	else
+  		y <- y + data$rate[[j]]/J
+    pop.total <- pop.total + data$pop[[j]]
+  }
+  rate.product <- demogdata(y, pop=pop.total, data$age, data$year, type=data$type, lambda=data$lambda,
+      label=data$label, name="product")
+  for (j in 1:J)
 	{
 		if(is.mortality)
 		{
@@ -48,15 +48,15 @@ coherentfdm <- function(data, order1=6, order2=6, ...)
 		rate.ratio[[j]]$rate[[1]][infrates] <- NA
 	}
 
-    # GM model
-    fdm.mean <- fdm(rate.product, series="product", order=order1, ...)
-    
-    # Ratio model
-    for (j in 1:J) 
-        fdm.ratio[[j]] <- fdm(rate.ratio[[j]], series=names(data$rate)[j], order=order2, ...)
-    names(fdm.ratio) <- names(data$rate)
-    
-    return(structure(list(product=fdm.mean, ratio=fdm.ratio), class="fdmpr"))
+  # GM model
+  fdm.mean <- fdm(rate.product, series="product", order=order1, ...)
+
+  # Ratio model
+  for (j in 1:J) 
+    fdm.ratio[[j]] <- fdm(rate.ratio[[j]], series=names(data$rate)[j], order=order2, ...)
+  names(fdm.ratio) <- names(data$rate)
+
+  return(structure(list(product=fdm.mean, ratio=fdm.ratio), class="fdmpr"))
 }
 
 
@@ -68,7 +68,7 @@ forecast.fdmpr <- function(object, h=50, level=80, K=100, drange=c(0.0,0.5), ...
   K <- min(K,ny)
 
   # GM model
-	fcast.mean <- forecast(object$product,method="arima",h=h,level=level,...)
+	fcast.mean <- forecast(object$product, method="arima", h=h, level=level, ...)
   # Make sure first coefficient is not I(1) with drift.
   #mod <- auto.arima(object$product$coeff[,2],d=2)
   #fcast.mean$coeff[[2]] <- forecast(mod, h=h, level=level, ...)

R/fdm.R

@@ -35,7 +35,7 @@ fdm <- function(data, series=names(data$rate)[1], order=6, ages=data$age, max.ag
 	nx <- length(data$age)
 	for(i in 1+(1:order))
 	{
-		if(sum(fit$basis[,i] > 0) < nx/2)
+		if(sum(na.omit(fit$basis[,i]) > 0) < nx/2)
 		{
 			fit$basis[,i] <- -fit$basis[,i]
 			fit$coeff[,i] <- -fit$coeff[,i]

R/hmd.R

@@ -1,10 +1,9 @@
 # Function to construct a mortality demogdata object from HMD
 hmd.mx <- function(country, username, password, label=country)
 {
-    require(RCurl)
     path <- paste("http://www.mortality.org/hmd/", country, "/STATS/", "Mx_1x1.txt", sep = "")
     userpwd <- paste(username, ":", password, sep = "")
-    txt <- getURL(path, userpwd = userpwd)
+    txt <- RCurl:::getURL(path, userpwd = userpwd)
     con <- textConnection(txt)
     mx <- try(read.table(con, skip = 2, header = TRUE, na.strings = "."),TRUE)
     close(con)
@@ -13,7 +12,7 @@ hmd.mx <- function(country, username, password, label=country)
 
     path <- paste("http://www.mortality.org/hmd/", country, "/STATS/", "Exposures_1x1.txt", sep = "")
     userpwd <- paste(username, ":", password, sep = "")
-    txt <- getURL(path, userpwd = userpwd)
+    txt <- RCurl:::getURL(path, userpwd = userpwd)
     con <- textConnection(txt)
     pop <- try(read.table(con, skip = 2, header = TRUE, na.strings = "."),TRUE)
     close(con)
@@ -46,12 +45,11 @@ hmd.mx <- function(country, username, password, label=country)
 }
 
 
-hmd.e0 = function(country, username, password)
+hmd.e0 <- function(country, username, password)
 {
-    require(RCurl)
     path <- paste("http://www.mortality.org/hmd/", country, "/STATS/", "E0per.txt", sep = "")
     userpwd <- paste(username, ":", password, sep = "")
-    txt <- getURL(path, userpwd = userpwd)
+    txt <- RCurl:::getURL(path, userpwd = userpwd)
     con <- textConnection(txt)
     lt <- try(read.table(con, skip = 2, header = TRUE, na.strings = "."), TRUE)
     close(con)

R/lca.R

@@ -140,10 +140,9 @@ lca <-  function(data,series=names(data$rate)[1],years=data$year, ages=data$age,
     {
         if(breakmethod=="bai")
         {
-            require(strucchange)
             x <- 1:m
             # Find breakpoints
-            bp <- breakpoints(kt ~ x)$breakpoints
+            bp <- strucchange:::breakpoints(kt ~ x)$breakpoints
             # Omit breakpoints less than minperiod from end
             bp <- bp[bp <= (m-minperiod)]
             bestbreak <- max(bp)