rocknRrr
Reputation: 421
State Space model code from Time Series analysis Springer book, and code error from optimization
I'm trying to apply State Space model code from the book, Time Series analysis and its applications from Springer book, to my data.
The code is as below:
num = length(jj)
A = cbind(1,1,0,0)
# Function to Calculate Likelihood
Linn =function(para){
Phi = diag(0,4)
Phi[1,1] = para[1]
Phi[2,]=c(0,-1,-1,-1)
Phi[3,]=c(0,1,0,0)
Phi[4,]=c(0,0,1,0)
cQ1 = para[2]; cQ2 = para[3] # sqrt q11 and q22
cQ = diag(0,4); cQ[1,1]=cQ1; cQ[2,2]=cQ2
cR = para[4]
kf = Kfilter0(num, jj, A, mu0, Sigma0, Phi, cQ, cR)
return(kf$like) }
# Initial Parameters
mu0 = c(.7,0,0,0); Sigma0 = diag(.04,4)
init.par = c(1.03,.1,.1,.5) # Phi[1,1], the 2 cQs and cR
# Estimation and Results
est = optim(init.par, Linn,NULL, method='BFGS', hessian=TRUE,
control=list(trace=1,REPORT=1))
SE = sqrt(diag(solve(est$hessian)))
u = cbind(estimate=est$par, SE)
rownames(u)=c('Phi11','sigw1','sigw2','sigv'); u
# Smooth
Phi = diag(0,4); Phi[1,1] = est$par[1]
Phi[2,]=c(0,-1,-1,-1); Phi[3,]=c(0,1,0,0); Phi[4,]=c(0,0,1,0)
cQ1 = est$par[2]; cQ2 = est$par[3]
cQ = diag(1,4); cQ[1,1]=cQ1; cQ[2,2]=cQ2
cR = est$par[4]
ks = Ksmooth0(num,jj,A,mu0,Sigma0,Phi,cQ,cR)
# Plots
Tsm = ts(ks$xs[1,,], start=1960, freq=4)
Ssm = ts(ks$xs[2,,], start=1960, freq=4)
p1 = 3*sqrt(ks$Ps[1,1,]); p2 = 3*sqrt(ks$Ps[2,2,])
par(mfrow=c(2,1))
plot(Tsm, main='Trend Component', ylab='Trend')
xx = c(time(jj), rev(time(jj)))
yy = c(Tsm-p1, rev(Tsm+p1))
polygon(xx, yy, border=NA, col=gray(.5, alpha = .3))
plot(jj, main='Data & Trend+Season', ylab='J&J QE/Share', ylim=c(-.5,17))
xx = c(time(jj), rev(time(jj)) )
yy = c((Tsm+Ssm)-(p1+p2), rev((Tsm+Ssm)+(p1+p2)) )
polygon(xx, yy, border=NA, col=gray(.5, alpha = .3))
# Forecast
n.ahead = 12
y = ts(append(jj, rep(0,n.ahead)), start=1960, freq=4) rmspe = rep(0,n.ahead)
x00 = ks$xf[,,num]; P00 = ks$Pf[,,num]
Q = t(cQ)%*%cQ
R = t(cR)%*%(cR)
for (m in 1:n.ahead){
xp = Phi%*%x00; Pp = Phi%*%P00%*%t(Phi)+Q
sig = A%*%Pp%*%t(A)+R; K = Pp%*%t(A)%*%(1/sig)
x00 = xp; P00 = Pp-K%*%A%*%Pp
y[num+m] = A%*%xp; rmspe[m] = sqrt(sig) }
plot(y, type='o', main='', ylab='J&J QE/Share', ylim=c(5,30),
xlim=c(1975,1984))
upp = ts(y[(num+1):(num+n.ahead)]+2*rmspe, start=1981, freq=4)
low = ts(y[(num+1):(num+n.ahead)]-2*rmspe, start=1981, freq=4)
xx = c(time(low), rev(time(upp)))
yy = c(low, rev(upp))
polygon(xx, yy, border=8, col=gray(.5, alpha = .3)) abline(v=1981, lty=3)
Then there's error for est (Estimation and Results).
> est = optim(init.par, Linn,NULL, method='BFGS', hessian=TRUE,
+ control=list(trace=1,REPORT=1))
Error in optim(init.par, Linn, NULL, method = "BFGS", hessian = TRUE, :
initial value in 'vmmin' is not finite
This is my data used for jj
jj<-c ( 102.95, 102.77, 103.06, 103.39, 102.98, 103.28, 102.64, 102.82, 102.14, 102.59, 103.12,
102.92, 103.65, NA, NA, 102.98, 102.86, 102.80, 102.30, 102.61, 102.83, 103.12,
102.78, 102.73, 102.46, 102.29, 102.73, 102.62, 102.56, 102.25, 102.22, 102.35, 102.45,
102.17, 102.05, 102.52, 101.87, NA, 102.16, 102.11, 102.14, 101.61, 101.96, 101.91,
102.08, 102.35, 101.95, 101.99, 101.70, 102.00, 101.93, 101.98, 101.94, 101.96, 101.77,
102.09, 102.22, 101.70, 102.03, 102.11, 101.95, 102.31, 101.91, 101.79, 103.18, 101.88,
102.02, NA, NA, 101.72, 102.06, 101.98, 102.10)
Why the error caused and anyone can run this code with the data I provide? The book doesn't provide data jj.
Upvotes: 1
Views: 100
Answers (1)
Tim-TU
Reputation: 408
the data jj is provided in the package astsa. If you have installed the package you can check it using astsa::jj. The class of astsa::jj is ts meaning that it is a time series object.
However, the error in your optimization is caused by the NA's.
I've attached a reprex where you can see that the example of the book seems to work.
And in a second example I've transformed your data into a ts object after eliminating the NA values. Optimization with optim worked for the adjusted data, but the algorithm has not reached convergence and stopped after 100 iterations. The computation of the standard errors throws a warning that NaNs were produced. The reason for that is that the variance-covariance matrix has a negative value on the diagonal (but this are the variances and they should be strictly positive).
Since I don't know the example taken from the book, you'll try and find out better starting values for your data.
In the attached reprex the forecast plot is at the end of the document (no idea why it is not in the right place), but that's not important.
library(astsa)
print(jj) # dataset inside package astsa
#> Qtr1 Qtr2 Qtr3 Qtr4
#> 1960 0.710000 0.630000 0.850000 0.440000
#> 1961 0.610000 0.690000 0.920000 0.550000
#> 1962 0.720000 0.770000 0.920000 0.600000
#> 1963 0.830000 0.800000 1.000000 0.770000
#> 1964 0.920000 1.000000 1.240000 1.000000
#> 1965 1.160000 1.300000 1.450000 1.250000
#> 1966 1.260000 1.380000 1.860000 1.560000
#> 1967 1.530000 1.590000 1.830000 1.860000
#> 1968 1.530000 2.070000 2.340000 2.250000
#> 1969 2.160000 2.430000 2.700000 2.250000
#> 1970 2.790000 3.420000 3.690000 3.600000
#> 1971 3.600000 4.320000 4.320000 4.050000
#> 1972 4.860000 5.040000 5.040000 4.410000
#> 1973 5.580000 5.850000 6.570000 5.310000
#> 1974 6.030000 6.390000 6.930000 5.850000
#> 1975 6.930000 7.740000 7.830000 6.120000
#> 1976 7.740000 8.910000 8.280000 6.840000
#> 1977 9.540000 10.260000 9.540000 8.729999
#> 1978 11.880000 12.060000 12.150000 8.910000
#> 1979 14.040000 12.960000 14.850000 9.990000
#> 1980 16.200000 14.670000 16.020000 11.610000
num = length(jj)
A = cbind(1,1,0,0)
# Function to Calculate Likelihood
Linn =function(para){
Phi = diag(0,4)
Phi[1,1] = para[1]
Phi[2,]=c(0,-1,-1,-1)
Phi[3,]=c(0,1,0,0)
Phi[4,]=c(0,0,1,0)
cQ1 = para[2]; cQ2 = para[3] # sqrt q11 and q22
cQ = diag(0,4); cQ[1,1]=cQ1; cQ[2,2]=cQ2
cR = para[4]
kf = Kfilter0(num, jj, A, mu0, Sigma0, Phi, cQ, cR)
return(kf$like) }
# Initial Parameters
mu0 = c(.7,0,0,0); Sigma0 = diag(.04,4)
init.par = c(1.03,.1,.1,.5) # Phi[1,1], the 2 cQs and cR
# Estimation and Results
est = optim(init.par, Linn,NULL, method='BFGS', hessian=TRUE,
control=list(trace=1, REPORT=1))
#> initial value 2.693644
#> iter 2 value -0.853526
#> iter 3 value -9.416505
#> iter 4 value -10.241752
#> iter 5 value -19.419809
#> iter 6 value -30.441188
#> iter 7 value -31.825543
#> iter 8 value -32.248413
#> iter 9 value -32.839918
#> iter 10 value -33.019870
#> iter 11 value -33.041749
#> iter 12 value -33.050583
#> iter 13 value -33.055492
#> iter 14 value -33.078152
#> iter 15 value -33.096870
#> iter 16 value -33.098405
#> iter 17 value -33.099018
#> iter 18 value -33.099385
#> iter 19 value -33.099498
#> iter 19 value -33.099498
#> final value -33.099498
#> converged
SE = sqrt(diag(solve(est$hessian)))
u = cbind(estimate=est$par, SE)
rownames(u)=c('Phi11','sigw1','sigw2','sigv'); u
#> estimate SE
#> Phi11 1.0350847657 0.00253645
#> sigw1 0.1397255477 0.02155155
#> sigw2 0.2208782663 0.02376430
#> sigv 0.0004655672 0.24174702
# Smooth
Phi = diag(0,4); Phi[1,1] = est$par[1]
Phi[2,]=c(0,-1,-1,-1); Phi[3,]=c(0,1,0,0); Phi[4,]=c(0,0,1,0)
cQ1 = est$par[2]; cQ2 = est$par[3]
cQ = diag(1,4); cQ[1,1]=cQ1; cQ[2,2]=cQ2
cR = est$par[4]
ks = Ksmooth0(num,jj,A,mu0,Sigma0,Phi,cQ,cR)
# Plots
Tsm = ts(ks$xs[1,,], start=1960, freq=4)
Ssm = ts(ks$xs[2,,], start=1960, freq=4)
p1 = 3*sqrt(ks$Ps[1,1,]); p2 = 3*sqrt(ks$Ps[2,2,])
par(mfrow=c(2,1))
plot(Tsm, main='Trend Component', ylab='Trend')
xx = c(time(jj), rev(time(jj)))
yy = c(Tsm-p1, rev(Tsm+p1))
polygon(xx, yy, border=NA, col=gray(.5, alpha = .3))
plot(jj, main='Data & Trend+Season', ylab='J&J QE/Share', ylim=c(-.5,17))
xx = c(time(jj), rev(time(jj)) )
yy = c((Tsm+Ssm)-(p1+p2), rev((Tsm+Ssm)+(p1+p2)) )
polygon(xx, yy, border=NA, col=gray(.5, alpha = .3))
# Forecast
n.ahead = 12
y = ts(append(jj, rep(0,n.ahead)), start=1960, freq=4)
rmspe = rep(0,n.ahead)
x00 = ks$xf[,,num]; P00 = ks$Pf[,,num]
Q = t(cQ)%*%cQ
R = t(cR)%*%(cR)
for (m in 1:n.ahead){
xp = Phi%*%x00; Pp = Phi%*%P00%*%t(Phi)+Q
sig = A%*%Pp%*%t(A)+R; K = Pp%*%t(A)%*%(1/sig)
x00 = xp; P00 = Pp-K%*%A%*%Pp
y[num+m] = A%*%xp; rmspe[m] = sqrt(sig) }
plot(y, type='o', main='', ylab='J&J QE/Share', ylim=c(5,30),
xlim=c(1975,1984))
upp = ts(y[(num+1):(num+n.ahead)]+2*rmspe, start=1981, freq=4)
low = ts(y[(num+1):(num+n.ahead)]-2*rmspe, start=1981, freq=4)
xx = c(time(low), rev(time(upp)))
yy = c(low, rev(upp))
polygon(xx, yy, border=8, col=gray(.5, alpha = .3))
abline(v=1981, lty=3)
# Example 2: New jj data:
jj_new <- c( 102.95, 102.77, 103.06, 103.39, 102.98, 103.28, 102.64, 102.82, 102.14, 102.59, 103.12,
102.92, 103.65, NA, NA, 102.98, 102.86, 102.80, 102.30, 102.61, 102.83, 103.12,
102.78, 102.73, 102.46, 102.29, 102.73, 102.62, 102.56, 102.25, 102.22, 102.35, 102.45,
102.17, 102.05, 102.52, 101.87, NA, 102.16, 102.11, 102.14, 101.61, 101.96, 101.91,
102.08, 102.35, 101.95, 101.99, 101.70, 102.00, 101.93, 101.98, 101.94, 101.96, 101.77,
102.09, 102.22, 101.70, 102.03, 102.11, 101.95, 102.31, 101.91, 101.79, 103.18, 101.88,
102.02, NA, NA, 101.72, 102.06, 101.98, 102.10)
jj <- ts(jj_new[!is.na(jj_new)])
num = length(jj)
A = cbind(1,1,0,0)
# Estimation and Results
est_new = optim(init.par, Linn,NULL, method='BFGS', hessian=TRUE,
control=list(trace=1, REPORT=1))
#> initial value 69066.649055
#> iter 2 value 13379.941416
#> iter 3 value 897.811170
#> iter 4 value 897.744837
#> iter 5 value 740.086083
#> iter 6 value 704.062900
#> iter 7 value 647.291087
#> iter 8 value 605.620617
#> iter 9 value 577.994261
#> iter 10 value 577.994023
#> iter 11 value 577.948157
#> iter 12 value 577.946077
#> iter 13 value 577.944346
#> iter 14 value 577.940792
#> iter 15 value 577.934726
#> iter 16 value 577.905521
#> iter 17 value 577.830290
#> iter 18 value 577.701675
#> iter 19 value 577.701480
#> iter 20 value 577.407053
#> iter 21 value 577.194089
#> iter 22 value 576.554158
#> iter 23 value 567.194439
#> iter 24 value 539.288477
#> iter 25 value 531.890806
#> iter 26 value 531.036336
#> iter 27 value 531.030643
#> iter 28 value 530.956134
#> iter 29 value 530.918352
#> iter 30 value 530.890127
#> iter 31 value 530.887646
#> iter 32 value 530.885713
#> iter 33 value 530.338301
#> iter 34 value 530.231457
#> iter 35 value 530.229078
#> iter 36 value 530.222573
#> iter 37 value 530.220160
#> iter 38 value 530.206317
#> iter 39 value 530.179131
#> iter 40 value 530.083081
#> iter 41 value 529.669557
#> iter 42 value 528.424740
#> iter 43 value 528.028025
#> iter 44 value 528.025772
#> iter 45 value 527.992265
#> iter 46 value 527.988454
#> iter 47 value 527.984442
#> iter 48 value 527.954276
#> iter 49 value 527.890590
#> iter 50 value 527.626120
#> iter 51 value 525.410701
#> iter 52 value 525.263544
#> iter 53 value 525.262081
#> iter 54 value 525.245961
#> iter 55 value 525.238839
#> iter 56 value 525.165379
#> iter 57 value 524.994817
#> iter 58 value 522.481418
#> iter 59 value 516.856037
#> iter 60 value 481.813924
#> iter 61 value 481.808421
#> iter 62 value 481.755703
#> iter 63 value 481.695417
#> iter 64 value 481.680106
#> iter 65 value 481.664492
#> iter 66 value 481.662879
#> iter 67 value 481.654499
#> iter 68 value 481.555816
#> iter 69 value 478.911532
#> iter 70 value 478.569138
#> iter 71 value 478.564163
#> iter 72 value 477.213754
#> iter 73 value 476.426603
#> iter 74 value 476.421185
#> iter 75 value 476.098605
#> iter 76 value 476.011245
#> iter 77 value 476.005941
#> iter 78 value 476.000627
#> iter 79 value 475.994552
#> iter 80 value 475.976022
#> iter 81 value 475.902453
#> iter 82 value 475.894178
#> iter 83 value 475.425270
#> iter 84 value 474.679311
#> iter 85 value 471.311778
#> iter 86 value 470.959303
#> iter 87 value 470.953483
#> iter 88 value 470.947660
#> iter 89 value 470.941837
#> iter 90 value 470.936012
#> iter 91 value 470.930187
#> iter 92 value 470.924361
#> iter 93 value 470.918534
#> iter 94 value 470.912707
#> iter 95 value 470.906878
#> iter 96 value 470.901049
#> iter 97 value 470.895219
#> iter 98 value 470.889388
#> iter 99 value 470.883557
#> iter 100 value 470.877724
#> final value 470.877724
#> stopped after 100 iterations
SE_new = sqrt(diag(solve(est_new$hessian))); SE_new
#> Warning in sqrt(diag(solve(est_new$hessian))): NaNs wurden erzeugt
#> [1] 25.22358 124.10763 NaN 103.69772
u_new = cbind(estimate=est_new$par, SE_new)
rownames(u_new)=c('Phi11','sigw1','sigw2','sigv'); u_new
#> estimate SE_new
#> Phi11 -0.8707774 25.22358
#> sigw1 -5.8028386 124.10763
#> sigw2 776.3957560 NaN
#> sigv 371.0135779 103.69772
Created on 2020-07-17 by the reprex package (v0.3.0)
Upvotes: 2
Related Questions
- Error in forecast and accuracy function in R
- Simulating ar(2) model generates an error message
- When ets() is used, why R is not responding and crashes?
- Question with using time series in R for forecasting via example
- Error when trying to apply forecasting methods on forecast object
- Implementing a Kalman Filter in MATLAB using 'ss'
- Using estimated state space parameters to forecast
- R Forecasting with covariates MARSS package
- R error on timeseries
- STLF function in the FORECAST package