search

helske / KFAS

KFAS: R Package for Exponential Family State Space Models
64 stars 17 forks source link

Difficulty with setting a1 and P1 when using SSMseasonal and reduced number of harmonics #77

Open bking124 opened 1 year ago

bking124 commented 1 year ago

When using SSMseasonal with trigonometric seasonality and a reduced number of harmonics, the current setup of the package is a bit counterintuitive if you want to set a1 and P1. As an example, I've provided code here for the sunspots.month data which is present in base R. This time series has an extended seasonal cycle, occurring about every 11 years. Hence if we want to model the seasonality, we will most likely need the Fourier representation, and likely don't need/want the full harmonics. So in the example below, I've written the KFAS code to limit to the first five harmonics.

library(KFAS)
plot(sunspot.month)
y <- as.vector(sunspot.month)
KFAS_mod <- SSModel(y ~ SSMtrend(degree=1, Q=list(diag(1)), a1=c(0), P1= c(diag(3))) +
                      SSMseasonal(period=12*11, Q=as.matrix(0), sea.type="trigonometric", harmonics=1:5))

This works well. However, if we want to set a1 and P1, we might think that we need to give 10 values, since there are 10 harmonics, i.e. we would write SSMseasonal(period=12*11, Q=as.matrix(0), sea.type="trigonometric", harmonics=1:5, a1=rep(0,10), P1= c(diag(3,10))

But this returns an error:

a1 must be a (m x 1) matrix where m is the number of states

To avoid the error, we have to provide an a1 and P1 that would cover the full set of harmonics, even though we are only using a reduced set. So, this works: SSMseasonal(period=12*11, Q=as.matrix(0), sea.type="trigonometric", harmonics=1:5, a1=rep(0,(12*11)-1), P1= c(diag(3,(12*11)-1))

As I said, I think this is a bit counterintuitive and in my mind makes more sense if it worked like the first example. This is definitely a corner case, but thought it worth raising nonetheless.

helske commented 1 year ago

Thanks, I added the harmonics argument at some point somewhat lazily, the SSMseasonal first constructs the system matrices using the full set of harmonics and then in the end subsets the matrices accordingly, so when checking a1 and P1 it the reduced dimensions are not accounted for. Indeed it would make sense that the checks worked in a way that allows your first example.