One step ahead predictions (dotted) and actual data (solid) bb=estBlackBox(my,g=3) #lag=?!!! tfplot(bb)
15 Series 1 0 5 10
1955
1960
1965
1970
1975
1980
Example: FKF
Example: FKF


#Shumway Example 6.12 newdat = read.table("f:/xzwu/2009/berkeley/ts/shumway/mydata/Newbold1.txt", comment.char="#") y = newdat[1:50,2] # quarterly inflation (column 2) z = newdat[1:50,3] # quarterly interest rates (column 3) plot(ts(y,start=c(1953,1),frequency=1)) points(ts(y,start=c(1953,1),frequency=1),type='p',pch=15) lines(ts(z,start=c(1953,1),frequency=1),lty=2) points(ts(z,start=c(1953,1),frequency=1),type='p',pch=18)
Step 2: Log Likelihood Function
• FKF: write function such as objective to create log likelihood function • Shumway: write function via Kfilter0, Kfilter1, Kfilter2 to create log likelihood function • dlm: write function such as buildFun
seriesNamesInput(my)="CPI" seriesNamesOutput(my)="BILLS“ bb=estBlackBox(my,g=3) #lag=?!!! tfplot(bb) bb$model #Giving F, G, H, K bb$estimates #giving $cov, $like, $pred attributes(bb)#model #看结果
• dse: skips this step
Step 3: Get estimates from result of step2
பைடு நூலகம்
• • • •
FKF: directly using R function optim Shumway directly using R function optim dlm: using function dlmMLE dse: using functions such as estVARXls, estVARXar, estSSfromVARX, estMaxLik, bft, estBlackBox, etc, directly
状态空间模型及例子
吴喜之
Models and Computation Steps for State-Space Model
Sources (for four packages): Package: FKF, related article (not for R, but having formulas): /p/dgr/kubcen/1998141.html (file3826.pdf) Package: dlm dlm.pdf (in R project or your R directory) Package: Shumway: /stoffer/tsa2/chap6.htm Package: dse In your R directory : dse-guide.pdf and article: http://www.bank-banque-canada.ca/pgilbert/gil93.pdf
ans <- fkf(a0 = c(0, 0), P0 = diag(c(10, 10)), dt = rbind(0, 0), ct = matrix(0), Tt = matrix(c(ar1, 1, ar2, 0), ncol = 2), Zt = cbind(1, 0), HHt = matrix(c(sigma^2, 0, 0, 0), ncol = 2), GGt = matrix(0), yt = rbind(a))
#library(dse) my1=TSdata(output=newdat[,c(3,2),drop = F]) my1=tframed(my1, list(start=c(1953,1), frequency=4)) bb1=estBlackBox(my1,g=3) tfplot(bb1) attributes(bb1)#model
Example: FKF
#Simulationdata ar1 <- 0.6 ar2 <- 0.2 ma1 <- -0.2 sigma <- sqrt(0.2) ## Sample from an ARMA(2, 1) process a <- arima.sim(model = list(ar = c(ar1, ar2), ma = ma1), n = n, innov = rnorm(n) * sigma)
– dlmModARMA , dlmModPoly, dlmModReg , dlmModSeas , dlmModTrig
• dse: using functions such as ARMA, SS (can transfer to each other) to build an empty model (but)
5
例子
ts(y, start = c(1953, 1), frequency = 1)
-2
-1
0
1
2
3
4
1960
1970 Time
1980
1990
2000
library(dse) my=TSdata(input= newdat[,3,drop = F],output=newdat[, 2,drop = F]) my=tframed(my, list(start=c(1953,1), frequency=4))
updating equations where
The MSE of
is Kalman gain
Kalman Smoother: with initial
The Lag-one Covariance Smoother
With initial condition
for t=n,n-1,…,2
Three “levels” of Shumway’s models
0 2
4
6
bb2$model #Giving F, H, K (one dimention) bb2$estimates #giving $cov, $like, $pred
8
10
12
14
1955
1960
1965
1970
1975
1980
bb4=estVARXls(my1,g=3)#two variables (TWO output) bb4$model #Giving A(L), B(L), C(L) bb4$estimates #giving $cov, $like, $pred
One step ahead predictions (dotted) and actual data (solid)
bb2=estBlackBox(my2,g=3) tfplot(bb2) attributes(bb2)#model
Series 1
bb3=estVARXls(my,g=1)#two variables (one input one output) bb3$model #Giving A(L), B(L), C(L) bb3$estimates #giving $cov, $like, $pred
Series 1
One step ahead predictions (dotted) and actual data (solid)
12
bb1$model #Giving F, H, K bb1$estimates #giving $cov, $like, $pred #看结果
0
4
8
1955
1960
Model (Shumway notation)
Transition equation and observation equation
• • • • • • •
xt : p×1 state vector N(m, S) yt : q×1 observation/measurement vector wt : p×1 system noise iid N(0,Q) vt : q×1 observation noise iid N(0,R) ut : r×1 fixed input At : q×p observation/measurement matrix F : p×p ; G : q×r, U: p×r, Q: p×p, R: q×q
library(dse) my=TSdata(input= newdat[,3,drop = F],output=newdat[, 2,drop = F]) my=tframed(my, list(start=c(1953,1), frequency=4)) seriesNamesInput(my)="CPI" seriesNamesOutput(my)="BILLS"