State space model with seasonal components

[Jennylai]

First thank you very much for solving my past questions.

I have a question of incoporating seasonal components into state space model and estimate them by DLM in RATS. My question is as follows:

I have monthly data of a series named FDI, and I select s=12, so if I want to incorporate a random walk component, a slope component and an AR(1) component into the state vector, I will have a state vector which has 14 elements, among these 14, 11 elements would be seasonal components.

The question I have is, if I only have three elements in the state vectors, to define matrices in 3 by 3 form is relatively easy, as my previous programme shows:

"environment noecho

*input datafile
open data FDIF_M_97_09.xls
data(format=xls,org=colulmns) / FDI
calendar(m) 1997:01
all 2009:03

*initialize the variables and matrices
nonlin v1 v2 v3
frml M3 = v1/sqrt(1+v1*v1)
frml Q1 = exp(v2)
frml Q3 = exp(v3)

dec rect z
dec rect h
*dec vect g
dec frml[rect] M
dec frml[symm] Q
dec frml[vect] f

compute z = ||1.0|0.0|1.0||
compute g = ||0.0|0.0|1.0||
compute h = 1.0

frml M = ||1.0,1.0,0.0|0.0,1.0,0.0|0.0,0.0,M3||
frml Q = ||Q1,0.0,0.0|0.0,0.0,0.0|0.0,0.0,Q3||
frml f = ||FDI||
"

And now, if I want to incorporate seasonal component into the state space model, I will have the matrices M, Q as 14 by 14, g as 14 by 1. Instead of writing them explicitly, which apparently will cost a huge space in programme file, I just wonder if there is an easy way to write these matrices in a compact way (as what we could have done in MATLAB or other software.), But I tried a lot and it does not work.

May I ask if you could give me some hint on this problem?

Thank you very much!

[TomDoan]

We have a procedure seasonaldlm.src which sets up the various matrices needed for the seasonal component of a state-space model. You can then combine those with the non-seasonal model that you've already built to create the overall SSM.

The following textbook examples use seasonaldlm:

Brockwell & Davis
itsfp281.prg

Durbin & Koopman
durkp162.prg
durkp167.prg
durkp230.prg
durkp233.prg

Tsay
tsayp535.prg

West & Harrison textbook
westp257.prg
westp318.prg
westp331.prg
westp387.prg
westp404.prg

[Jennylai]

Dear Tom ,

Thank you very much! I will go and check it out.

[Jennylai]

Dear Tom,

May I ask about this "phis" in the durkp162.r file in RATS?

there is a line in durkp162.r: compute phis=-4.0.

May I ask how you come up with this idea of phis = -4.0?

Thank you very much!

[Jennylai]

Dear Tom,

I have another question. If I want to expand the program file localdlm.r, to include the AR(1) component. But there is an unknown parameter related to the A matrix if AR(1) component is being included. In this case, is it still possible to 'glue' the Al matrix with As matrix from seasonal components?

[Jennylai]

Dear Tom,

I tried to incorporate AR(1) part into the state vectors, and I just revised a little bit of the program durkp162 by adding M3 into Nonlin:
"
nonlin phil phis v1
frml m3=v1/sqrt(1+v1*v1)
.
.
.
dlm(method=bfgs,start=(a=al~\M3~\as,sw=exp(2*phil)*swl~\exp(2*phis)*sws),$
y=FDI,a=a,c=c,exact,scale,sv=1.0,sw=sw,type=filter) / xfilter

"

But it does not work. May I ask how to incorporate this AR(1) part into the model? Thank you!

[TomDoan]

Dear Tom,

May I ask about this "phis" in the durkp162.r file in RATS?

there is a line in durkp162.r: compute phis=-4.0.

May I ask how you come up with this idea of phis = -4.0?

Thank you very much!

Probably some experimentation. It needs to be "small" relative to the other component standard deviations. (It's the log of the ratio of the seasonal shock s.d. to the observation error s.d.).

[TomDoan]

Dear Tom,

I tried to incorporate AR(1) part into the state vectors, and I just revised a little bit of the program durkp162 by adding M3 into Nonlin:

nonlin phil phis v1
frml m3=v1/sqrt(1+v1*v1)
.
.
.
dlm(method=bfgs,start=(a=al~\M3~\as,sw=exp(2*phil)*swl~\exp(2*phis)*sws),$
y=FDI,a=a,c=c,exact,scale,sv=1.0,sw=sw,type=filter) / xfilter

"

But it does not work. May I ask how to incorporate this AR(1) part into the model? Thank you!

Code: Select all

frml m3=v1/sqrt(1+v1*v1)

You need a space before the equal sign.

Code: Select all

dlm(method=bfgs,start=(a=al~\M3~\as,sw=exp(2*phil)*swl~\exp(2*phis)*sws),$
   y=FDI,a=a,c=c,exact,scale,sv=1.0,sw=sw,type=filter) / xfilter

The dimension of SW is going to be wrong here. swl has the shocks to the linear trend model, sws has the shocks to the seasonal, but you need the variance to the shock to the AR in between them.

Also, make sure that your "C" matrix adds in the AR component.

[Jennylai]

Dear Tom,

Thank you very much for your kind reply. It really helps a lot!

By the way, I have another question regarding the LOCALDLM.src file.

In this procedure, under the option "type=trend", you wrote a line " if %defined(sw) /n compute sw=||0.0|0.0,1.0||"

I am a little bit puzzled at this line. I thought it might be that " sw =||1.0|0.0,0.0||", because the level(random walk component) has an error term that has an variance, and we sometimes suppress the error term in the slope(usually called beta) to be zero and the variance will accordingly be zero.

What do you think?

[TomDoan]

Dear Tom,

Thank you very much for your kind reply. It really helps a lot!

By the way, I have another question regarding the LOCALDLM.src file.

In this procedure, under the option "type=trend", you wrote a line " if %defined(sw) /n compute sw=||0.0|0.0,1.0||"

I am a little bit puzzled at this line. I thought it might be that " sw =||1.0|0.0,0.0||", because the level(random walk component) has an error term that has an variance, and we sometimes suppress the error term in the slope(usually called beta) to be zero and the variance will accordingly be zero.

What do you think?

Although it's possible to include independent shocks in both, it's rather common for the variance in the level to come in zero - the combination of three independent shocks (trend rate, level and measurement error) is too many to identify separately.

[Jennylai]

Dear Tom,

Thank you !

I have another question concerning the command in durkp162.r
" compute phis=-4.0 "

My question is this: If the variance of the irregular part in the measurement equation is being concentrated out of the likelihood function, then the variance of the seasonal error term is porportional to the variance of the irregular part. If this is true, then shall we write " compute phis=log(sqrt(vs/vi))" instead of the " compute phis=-4.0", vs stands for the variance for the error term of the seasonal components.

[Jennylai]

Dear Tom,

Just forget about my last post, I worked it out myself~~

Thank you very much for your kind reply all the time!

[Jennylai]

Dear Tom, About the AR(1) component, you suggested last time that:

Dear Tom,

I tried to incorporate AR(1) part into the state vectors, and I just revised a little bit of the program durkp162 by adding M3 into Nonlin:

nonlin phil phis v1
frml m3=v1/sqrt(1+v1*v1)
.
.
.
dlm(method=bfgs,start=(a=al~\M3~\as,sw=exp(2*phil)*swl~\exp(2*phis)*sws),$
y=FDI,a=a,c=c,exact,scale,sv=1.0,sw=sw,type=filter) / xfilter

"

But it does not work. May I ask how to incorporate this AR(1) part into the model? Thank you!

Code: Select all

frml m3=v1/sqrt(1+v1*v1)

You need a space before the equal sign.

Code: Select all

dlm(method=bfgs,start=(a=al~\M3~\as,sw=exp(2*phil)*swl~\exp(2*phis)*sws),$
   y=FDI,a=a,c=c,exact,scale,sv=1.0,sw=sw,type=filter) / xfilter

The dimension of SW is going to be wrong here. swl has the shocks to the linear trend model, sws has the shocks to the seasonal, but you need the variance to the shock to the AR in between them.

Also, make sure that your "C" matrix adds in the AR component.

Dear Tom,

I have another question concerning the AR(1) component being incorporated into the State Space Model.

I've applied what you mentioned last time to the programme, and the revised programme looks like this:

open plot originalFDI.rgf
graph(footer="Figure 1 Monthly FDI inflows to China")
# FDI
close plot
*
* Generate the basic seasonal and local level component models
* SeasonalDLM and LocalDLM create the "A", "C" and "SW" matrices (in RATS notation)
*
@SeasonalDLM(type=fourier) as cs sws
@LocalDLM(type=trend) al cl swl
*
* Add the AR(1) part in to the matrix "C"
compute ca=||1.0||
*
* 'Glue' cl,ca and cs together. As Matrix "A" and "SW" will have unknown parameters,
* we will glue them inside the 'DLM' command.
*
compute c=cl~~ca~~cs
*
* Generate the 'SW' matrix for AR(1) component.
compute swa=||1.0||
*
*
nonlin phil phis rov phia
*
* The success of the estimation is quite sensitive to the initial guess values.
* The following is a bit cruder than the guess values computed using the
* reference in the book, but should work in practice
*
* (a) get rid of the seasonal by regression. (This corresponds to
* a zero variance for that component)
*
filter(remove=seasonal) FDI / deseas
*
* (b) do a short span filter on the deseasonalized data to get an
* estimate of the local trend
*
filter(type=centered,width=5) deseas / trend
*
* (c) compute the approximate irregular component and its variance
*
set irreg = deseas-trend
stats(noprint) irreg
compute vi=%variance
*
* (d) compute the first difference of the local level component and
* its variance
*
set difft = trend-trend{1}
stats(noprint) difft
compute vl=%variance
*
* Because we're concentrating out the irregular variance, the "phil" parameter
* needs to be initialized based upon the ratio. We'll make "phis" small.
*
* Set the initial values for 'phil','phis',and 'v3','va'.
*
compute v3=0
compute va=0
*
compute phil=log(sqrt(vl/vi))
compute phis=-4.0
compute rov=v3/sqrt(1+v3*v3)
compute phia=log(sqrt(va/vi))
*
* Generage the 'G' option for 'Exact' initialization to perform.
*
compute g=%zeros(1,%rows(c))
compute g(1,3)=1.0
*
* Because the A and SW matrix for the complete model depend upon the parameters being
* estimated, we couldn't "glue" it together earlier.
*
dlm(method=bfgs,start=(a=al~~rov~~as,sw=exp(2*phil)*swl~\exp(2*phia)*swa~\exp(2*phis)*sws),$
y=FDI,a=a,c=c,exact,g=g,scale,sv=1.0,sw=sw,type=filter) / xfilter pfilter

However, it does not work in RATS, and RATS reports the error like this:

dlm(method=bfgs,start=(a=al~~rov~~as,<<<<
## SX20. Expected ) Here
>>>>sw=exp(2*phil)*swl~\exp(2*phia)*swa~\exp(2*phis)*sws),$
y=FDI,a=a,c=c,exact,g=g,scale,sv=1.0,sw=sw,type=filter) / xfilter pfilter

Is it the case that RATS does not allow two expression inside the 'START' option in 'DLM' command?

[TomDoan]

However, it does not work in RATS, and RATS reports the error like this:

dlm(method=bfgs,start=(a=al~~rov~~as,<<<<
## SX20. Expected ) Here
>>>>sw=exp(2*phil)*swl~\exp(2*phia)*swa~\exp(2*phis)*sws),$
y=FDI,a=a,c=c,exact,g=g,scale,sv=1.0,sw=sw,type=filter) / xfilter pfilter
Is it the case that RATS does not allow two expression inside the 'START' option in 'DLM' command?

If you put a % in front of the (...,...) expression, it will make it possible to do that. That is:

Code: Select all

dlm(method=bfgs,start=%(a=al~~rov~~as,sw=exp(2*phil)*swl~\exp(2*phia)*swa~\exp(2*phis)*sws),$
y=FDI,a=a,c=c,exact,g=g,scale,sv=1.0,sw=sw,type=filter) / xfilter pfilter

The alternative (which is probably more readable) is to create a short function that does the calculation:

Code: Select all

function MyStart
compute a=al~~rov~~as
compute sw=sw=exp(2*phil)*swl~\exp(2*phia)*swa~\exp(2*phis)*sws
end MyStart
dlm(method=bfgs,start=MyStart(),$
y=FDI,a=a,c=c,exact,g=g,scale,sv=1.0,sw=sw,type=filter) / xfilter pfilter

[Jennylai]

Dear Tom,

Thank you very much for your help all the time!