State Space Time Series analysis Springer book 中的模型代码，以及优化中的代码错误
State Space model code from Time Series analysis Springer book, and code error from optimization
我正在尝试将 State Space 书中的模型代码、Springer 书中的时间序列分析及其应用程序应用到我的数据中。
代码如下：
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)
那么 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
这是我用于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)
为什么会导致错误，任何人都可以 运行 使用我提供的数据编写此代码？本书没有提供数据jj.
数据jj在包astsa中提供。如果您已经安装了该软件包，您可以使用 astsa::jj 检查它。 astsa::jj的class是ts表示是时间序列对象
但是，您的优化错误是由 NA 引起的。
我附上了一个 reprex，您可以在其中看到书中的示例似乎有效。
在第二个示例中，我在消除 NA 值后将您的数据转换为 ts 对象。使用optim优化对调整后的数据有效，但算法未达到收敛并在100次迭代后停止。标准错误的计算会发出警告，指出已生成 NaNs。原因是方差-协方差矩阵在对角线上有一个负值（但这是方差，它们应该严格为正）。
由于我不知道书中的示例，您将尝试为您的数据找到更好的起始值。
在所附的 reprex 中，预测图位于文档的末尾（不知道为什么它不在正确的位置），但这并不重要。
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
由 reprex package (v0.3.0)
于 2020-07-17 创建
我正在尝试将 State Space 书中的模型代码、Springer 书中的时间序列分析及其应用程序应用到我的数据中。
代码如下：
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)
那么 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
这是我用于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)
为什么会导致错误，任何人都可以 运行 使用我提供的数据编写此代码？本书没有提供数据jj.
数据jj在包astsa中提供。如果您已经安装了该软件包，您可以使用 astsa::jj 检查它。 astsa::jj的class是ts表示是时间序列对象
但是，您的优化错误是由 NA 引起的。
我附上了一个 reprex，您可以在其中看到书中的示例似乎有效。
在第二个示例中，我在消除 NA 值后将您的数据转换为 ts 对象。使用optim优化对调整后的数据有效，但算法未达到收敛并在100次迭代后停止。标准错误的计算会发出警告，指出已生成 NaNs。原因是方差-协方差矩阵在对角线上有一个负值（但这是方差，它们应该严格为正）。
由于我不知道书中的示例，您将尝试为您的数据找到更好的起始值。
在所附的 reprex 中，预测图位于文档的末尾（不知道为什么它不在正确的位置），但这并不重要。
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
由 reprex package (v0.3.0)
于 2020-07-17 创建