1
$\begingroup$

I think I found the solution myself but would need some verification by an expert. To see my solution you can skip the start and switch to the end of my question. My problem is now: How do I get a h-step forecast from there?

Dear Cross Validated Community,

I am trying to reproduce a paper in which a state space model of the form,

$p_t=\beta_ty_t$

$\beta_t=\alpha + \gamma\beta_{t-1} + \epsilon_{\beta, t}$

is used to "rescale" $y$ to $p$ ($y$ comes from another estimation, it retraces $p$ but is not on the same scale; unfortunately I am completely new to time series techniques). I tried it with the dlm package in R package, but I read in another post on Cross Validated that only the FKF package in R allows to explicitly define intercepts in the transition equation. So, although this package is unfortunately not very well documented, I tried my best. Here is my code:

library(xts)
library(zoo)
library(FKF)

rescp = structure(c(-1.48383556632677, -1.37794286288647, -0.223968294200336, 
1.91458717447093, 2.81412312320836, 2.62111138897671,     2.15133635786653, 
2.28617143932799, 2.20677075253045, 2.00354892921858,     1.12156059406863, 
1.3407784804821, 1.15017753951125, 1.11831224723318, 1.16669514893146, 
1.08454477660035, 1.43779302221793, 1.68285146766718, 2.12765957446808, 
2.61832134822156, 3.06205526142694, 3.48796539187262,    3.51843240530573, 
3.59091013533697, 3.74398296646063, 3.84881456877791,      3.55011756112037, 
3.44596748485815, 3.02668819623354, 2.98940076328974,      2.87556781310837, 
2.61024730546439, 2.28153173257163, 1.74145836020083,      1.68341183325105, 
1.40405200981283, 1.69363328349891, 1.95344897563043,      2.17361591809921, 
1.7807067358699, 1.76979462457901, 1.6179382593004, 1.97910506603588, 
1.47640977690002, 1.09088370307018, 1.40087855410102,      1.73077427764121, 
1.92314422706751, 1.52319810709362, 1.15489130434782,      -0.424786746371165, 
-0.405158849961268, -0.28598901718276, -0.325954647412152,      -0.434872061884164, 
-0.29270970418345, -0.294099020293278, -0.300697596065055,  -0.207938375799877, 
-0.172517121486599, -0.276722670510142, -0.289548544267883, -0.372364783542654, 
-0.369099557169152, -0.25999808351562, -0.233080646843824, -0.248548002943304, 
-0.0626764955490915, -0.00789927715314176, 0.196226273709107, 
0.207957770647136, 0.153011835302311, -0.0195874547653262,      -0.109097168921427, 
-0.088909879547551, -0.130316061633812, -0.122895889350873,      0.0103480904803911, 
-0.0548531386506467, 0.0452201814900919, 0.0148692080447909, 
0.0205034263264671, -0.0237959158194401, -0.0706051207969854, 
-0.207429686705238, -0.0672502087028582, 0.0246514539902095,   
-0.052667262253806, 0.0356379822532804, -0.0641823305528641, 
-0.0744098751098921, 0.0134397296338576, 0.0691642036631347, 
0.0967877697374845, -0.0193097496328191, -0.100499719682672, 
-0.00747754124178514, 0.0120458898383589, 0.130091772776046, 
0.0224959290304737), class = c("xts", "zoo"), .indexCLASS =         c("POSIXct", 
"POSIXt"), .indexTZ = "", tclass = c("POSIXct", "POSIXt"), tzone = "",     index = structure(c(1249077600, 
1251756000, 1254348000, 1257030000, 1259622000, 1262300400,      1264978800, 
1267398000, 1270072800, 1272664800, 1275343200, 1277935200,     1280613600, 
1283292000, 1285884000, 1288566000, 1291158000, 1293836400,     1296514800, 
1298934000, 1301608800, 1304200800, 1306879200, 1309471200,     1312149600, 
1314828000, 1317420000, 1320102000, 1322694000, 1325372400,      1328050800, 
1330556400, 1333231200, 1335823200, 1338501600, 1341093600,      1343772000, 
1346450400, 1349042400, 1351724400, 1354316400, 1356994800,      1359673200, 
1362092400, 1364767200, 1367359200, 1370037600, 1372629600,       1375308000, 
1377986400), tzone = "", tclass = c("POSIXct", "POSIXt")), .Dim =      c(50L, 
2L), .Dimnames = list(NULL, c("cp", "mu")))



#measurement equation: yt = ct + Zt at + Gt et 
#transition equation:  at+1 = dt + Tt at + Ht mut  

#dt either a m × n (time-varying) or a m × 1 (constant) matrix.
#Tt either a m × m × n or a m × m × 1 array.
#HHt either a m × m × n or a m × m × 1 array.
#ct either a d × n or a d × 1 matrix.
#Zt either a d × m × n or a d × m × 1 array.
#GGt either a d × d × n or a d × d × 1 array.
#yt a d × n matrix.

n   = nrow(rescp)
d   = 1
m   = 1
yt  = matrix(rescp[,"cp"],d,n)
ct  = matrix(0)
Zt  = array(rescp[,"mu"],c(d,m,n))
GGt = array(0,c(d,d,1))

kaf <- function(par,...){-fkf(Tt=matrix(par[1]),dt=matrix(par[2]),HHt=array(par[3],c(m,m,n)),...,check.input = T)$logLik}

par <- c(0.1,0.1,0.1)
fit <- optim(par,kaf,yt=yt,Zt=Zt,ct=ct,GGt=GGt,a0=a0,P0=P0,hessian=TRUE)
res <- fkf(Tt=matrix(fit$par[1]),dt=matrix(fit$par[2]),HHt=array(fit$par[3],c(m,m,n)),yt=yt,Zt=Zt,ct=ct,GGt=GGt,a0=a0,P0=P0,check.input = T)

The model runs and converges, but it would be very nice if someone could help me, as I am almost certain that I mis-specified something (the results are very strange). The code is supposed to run if you just copy and paste it into R. Departing from this model I would need a 12 month forecast at each point of time for $y$ (the data is monthly). How could I achieve that? I am grateful for every comment!

EDIT: Likely solution

I think I solved the problem myself. Below is the code I used. I still would like to have verification from an expert though. And I have an additional question: How do I get a h-step forecast from there?

n   = nrow(rescp)
d   = 1
m   = 1
yt  = matrix(rescp[,"cp"],d,n)
ct  = matrix(0)
Zt  = array(rescp[,"mu"],c(d,m,n))
GGt = array(0,c(d,d,1))
a0  = yt[,1]
P0  = matrix(10)

kaf <- function(par,...){-fkf(Tt=matrix(par[1]),dt=matrix(par[2]),HHt=array(par[3],c(m,m,n)),...,check.input = T)$logLik}
    par <- c(0.1,0.1,0.1)
    fit <- optim(par,kaf,yt=yt,Zt=Zt,ct=ct,GGt=GGt,a0=a0,P0=P0,hessian=TRUE)
    res <- fkf(Tt=matrix(fit$par[1]),dt=matrix(fit$par[2]),HHt=array(fit$par[3],c(m,m,n)),yt=yt,Zt=Zt,ct=ct,GGt=GGt,a0=a0,P0=P0,check.input = T)
$\endgroup$
  • $\begingroup$ Can you give the reference of the paper? $\endgroup$ – javlacalle Jun 5 '15 at 15:10
  • $\begingroup$ Shouldn't par[3] be related to the variance of the transition equation, HHt instead of GGt. Similarly, I think that ct should be zero and dt be defined as the intercept to be estimated, $\alpha$. $\endgroup$ – javlacalle Jun 5 '15 at 15:11
  • $\begingroup$ @javlacalle Thank your for your reply! Sure, the paper is: papers.ssrn.com/sol3/papers.cfm?abstract_id=2359500 (page 23, first paragraph; unfortunately I only found the working paper version for free). $\endgroup$ – chameau13 Jun 5 '15 at 15:16
  • $\begingroup$ @javlacalle Regarding your first comment. I just saw these mistakes myself. I will edit the question within the next twenty minutes. $\endgroup$ – chameau13 Jun 5 '15 at 15:19
  • $\begingroup$ @javlacalle Sorry to abuse this a some sort of chat, but would you please be so nice to look at my question again? I added a solution which I found myself and I think it is pretty much correct. However I do I get a h-step forecast for y from there? It would be very nice if you could help me! $\endgroup$ – chameau13 Jun 7 '15 at 21:52

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Browse other questions tagged or ask your own question.