Forecast accuracy in R

I have followed the instructions in this StMoMo package document to fit Lee Carter to mortality data for Canada.

The next step in my project is to measure the forecast accuracy of the Lee Carter model when fit to this Canadian data.

To do this, I tried to use accuracy() but ran into an error since my Lee Carter fit is of class "fitStMoMo" and not of class "forecast" or a time series.

Is there an alternative forecast accuracy function that I can use on "fitStMoMo" objects that will calculate Mean Error, Root Mean Squared Error, Mean Absolute Error, Mean Percentage Error, Mean Absolute Percentage Error and Mean Absolute Scaled Error for me?

Reprex: Reprex created using EWMaleData as used in the StMoMo document to flag the error specifically:

keep.source.pkgs = TRUE
install.packages("StMoMo")
#> 
#> The downloaded binary packages are in
#>  /var/folders/0k/cmrjt5fj0993lz39d1b8rmfh0000gp/T//RtmpvXMphw/downloaded_packages
install.packages("demography")
#> 
#> The downloaded binary packages are in
#>  /var/folders/0k/cmrjt5fj0993lz39d1b8rmfh0000gp/T//RtmpvXMphw/downloaded_packages
install.packages("forecast")
#> 
#> The downloaded binary packages are in
#>  /var/folders/0k/cmrjt5fj0993lz39d1b8rmfh0000gp/T//RtmpvXMphw/downloaded_packages
library("StMoMo")
#> Loading required package: gnm
#> Loading required package: forecast
library("demography")
#> This is demography 1.20
library("forecast")
constLC <- function(ax, bx, kt, b0x, gc, wxt, ages){
+     +     + c1 <- mean(kt[1, ], na.rm = TRUE)
+     +     + c2 <- sum(bx[, 1], na.rm = TRUE)
+     +     + list(ax = ax + c1 * bx, bx = bx / c2, kt = c2 * (kt - c1))
}
LC <- StMoMo(link = "logit", staticAgeFun = TRUE, periodAgeFun = "NP",
constFun = constLC)
LC <- lc(link = "logit")
LC$gnmFormula
#> [1] "D/E ~ -1 + offset(o) + factor(x) + Mult(factor(x), factor(t), inst = 1)"
EWMaleData
#> Mortality data for England and Wales
#>     Series:  male
#>     Years: 1961 - 2011
#>     Ages:  0 - 100
#>     Exposure:  central
EWMaleIniData <- central2initial(EWMaleData)
ages.fit <- 55:89
wxt <- genWeightMat(ages = ages.fit, years = EWMaleIniData$years,
clip = 3)
LCfit <- fit(LC, data = EWMaleIniData, ages.fit = ages.fit, wxt = wxt)
#> StMoMo: The following cohorts have been zero weigthed: 1872 1873 1874 1954 1955 1956 
#> StMoMo: Start fitting with gnm
#> Initialising
#> Running start-up iterations..
#> Running main iterations.....
#> Done
#> StMoMo: Finish fitting with gnm
LCfor <- forecast(LCfit, h = 50)
class(LCfit)
#> [1] "fitStMoMo"
class(LCfor)
#> [1] "forStMoMo"
accuracy(LCfit)
#> Error in accuracy.default(LCfit): First argument should be a forecast object or a time series.
accuracy(LCfor)
#> Error in accuracy.default(LCfor): First argument should be a forecast object or a time series.
Created on 2018-11-21 by the reprex package (v0.2.1)

Take a look at

help(forecast

and try it without the h=50 argument

> LCfor <- forecast(LCfit)
> LCfor
Stochastic Mortality Model forecast
Call: forecast.fitStMoMo(object = LCfit)

Binomial model with predictor: logit q[x,t] = a[x] + b1[x] k1[t]

kt model: mrwd
Jump-off method: fit
Data:  England and Wales
Series:  male
Years in forecast: 2012 - 2061
Ages in forecast: 55 - 89 

h is the number of periods in the forecast and it's set by default by formula. If you have to have 50, rather than 49, you'll need some thinkering somewhere in the arguments

1 Like

This topic was automatically closed 21 days after the last reply. New replies are no longer allowed.