Send output to:
Browser Blue - Charts White
Browser Black/White
CSV
Data:
6155 6490 6285 6450 6240 6375 6100 5600 5505 5155 4720 4645 5210 5580 5830 6395 7265 7460 6790 6810 7040 7385 7170 6620 6485 6035 6165 6320 6430 6785 6895 8030 8115 8915 9630 9250 9535 8320 7070 6375 6535 6480 6925 7020
Seasonal period
12
12
1
2
3
4
5
6
7
8
9
10
11
12
Number of Forecasts
12
12
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
Algorithm
BFGS
BFGS
L-BFGS-B
Chart options
R Code
require('stsm') require('stsm.class') require('KFKSDS') par1 <- as.numeric(par1) par2 <- as.numeric(par2) nx <- length(x) x <- ts(x,frequency=par1) m <- StructTS(x,type='BSM') print(m$coef) print(m$fitted) print(m$resid) mylevel <- as.numeric(m$fitted[,'level']) myslope <- as.numeric(m$fitted[,'slope']) myseas <- as.numeric(m$fitted[,'sea']) myresid <- as.numeric(m$resid) myfit <- mylevel+myseas mm <- stsm.model(model = 'BSM', y = x, transPars = 'StructTS') fit2 <- stsmFit(mm, stsm.method = 'maxlik.td.optim', method = par3, KF.args = list(P0cov = TRUE)) (fit2.comps <- tsSmooth(fit2, P0cov = FALSE)$states) m2 <- set.pars(mm, pmax(fit2$par, .Machine$double.eps)) (ss <- char2numeric(m2)) (pred <- predict(ss, x, n.ahead = par2)) mylagmax <- nx/2 bitmap(file='test2.png') op <- par(mfrow = c(2,2)) acf(as.numeric(x),lag.max = mylagmax,main='Observed') acf(mylevel,na.action=na.pass,lag.max = mylagmax,main='Level') acf(myseas,na.action=na.pass,lag.max = mylagmax,main='Seasonal') acf(myresid,na.action=na.pass,lag.max = mylagmax,main='Standardized Residals') par(op) dev.off() bitmap(file='test3.png') op <- par(mfrow = c(2,2)) spectrum(as.numeric(x),main='Observed') spectrum(mylevel,main='Level') spectrum(myseas,main='Seasonal') spectrum(myresid,main='Standardized Residals') par(op) dev.off() bitmap(file='test4.png') op <- par(mfrow = c(2,2)) cpgram(as.numeric(x),main='Observed') cpgram(mylevel,main='Level') cpgram(myseas,main='Seasonal') cpgram(myresid,main='Standardized Residals') par(op) dev.off() bitmap(file='test1.png') plot(as.numeric(m$resid),main='Standardized Residuals',ylab='Residuals',xlab='time',type='b') grid() dev.off() bitmap(file='test5.png') op <- par(mfrow = c(2,2)) hist(m$resid,main='Residual Histogram') plot(density(m$resid),main='Residual Kernel Density') qqnorm(m$resid,main='Residual Normal QQ Plot') qqline(m$resid) plot(m$resid^2, myfit^2,main='Sq.Resid vs. Sq.Fit',xlab='Squared residuals',ylab='Squared Fit') par(op) dev.off() bitmap(file='test6.png') par(mfrow = c(3,1), mar = c(3,3,3,3)) plot(cbind(x, pred$pred), type = 'n', plot.type = 'single', ylab = '') lines(x) polygon(c(time(pred$pred), rev(time(pred$pred))), c(pred$pred + 2 * pred$se, rev(pred$pred)), col = 'gray85', border = NA) polygon(c(time(pred$pred), rev(time(pred$pred))), c(pred$pred - 2 * pred$se, rev(pred$pred)), col = ' gray85', border = NA) lines(pred$pred, col = 'blue', lwd = 1.5) mtext(text = 'forecasts of the observed series', side = 3, adj = 0) plot(cbind(x, pred$a[,1]), type = 'n', plot.type = 'single', ylab = '') lines(x) polygon(c(time(pred$a[,1]), rev(time(pred$a[,1]))), c(pred$a[,1] + 2 * sqrt(pred$P[,1]), rev(pred$a[,1])), col = 'gray85', border = NA) polygon(c(time(pred$a[,1]), rev(time(pred$a[,1]))), c(pred$a[,1] - 2 * sqrt(pred$P[,1]), rev(pred$a[,1])), col = ' gray85', border = NA) lines(pred$a[,1], col = 'blue', lwd = 1.5) mtext(text = 'forecasts of the level component', side = 3, adj = 0) plot(cbind(fit2.comps[,3], pred$a[,3]), type = 'n', plot.type = 'single', ylab = '') lines(fit2.comps[,3]) polygon(c(time(pred$a[,3]), rev(time(pred$a[,3]))), c(pred$a[,3] + 2 * sqrt(pred$P[,3]), rev(pred$a[,3])), col = 'gray85', border = NA) polygon(c(time(pred$a[,3]), rev(time(pred$a[,3]))), c(pred$a[,3] - 2 * sqrt(pred$P[,3]), rev(pred$a[,3])), col = ' gray85', border = NA) lines(pred$a[,3], col = 'blue', lwd = 1.5) mtext(text = 'forecasts of the seasonal component', side = 3, adj = 0) dev.off() load(file='createtable') a<-table.start() a<-table.row.start(a) a<-table.element(a,'Structural Time Series Model -- Interpolation',6,TRUE) a<-table.row.end(a) a<-table.row.start(a) a<-table.element(a,'t',header=TRUE) a<-table.element(a,'Observed',header=TRUE) a<-table.element(a,'Level',header=TRUE) a<-table.element(a,'Slope',header=TRUE) a<-table.element(a,'Seasonal',header=TRUE) a<-table.element(a,'Stand. Residuals',header=TRUE) a<-table.row.end(a) for (i in 1:nx) { a<-table.row.start(a) a<-table.element(a,i,header=TRUE) a<-table.element(a,x[i]) a<-table.element(a,mylevel[i]) a<-table.element(a,myslope[i]) a<-table.element(a,myseas[i]) a<-table.element(a,myresid[i]) a<-table.row.end(a) } a<-table.end(a) table.save(a,file='mytable.tab') a<-table.start() a<-table.row.start(a) a<-table.element(a,'Structural Time Series Model -- Extrapolation',4,TRUE) a<-table.row.end(a) a<-table.row.start(a) a<-table.element(a,'t',header=TRUE) a<-table.element(a,'Observed',header=TRUE) a<-table.element(a,'Level',header=TRUE) a<-table.element(a,'Seasonal',header=TRUE) a<-table.row.end(a) for (i in 1:par2) { a<-table.row.start(a) a<-table.element(a,i,header=TRUE) a<-table.element(a,pred$pred[i]) a<-table.element(a,pred$a[i,1]) a<-table.element(a,pred$a[i,3]) a<-table.row.end(a) } a<-table.end(a) table.save(a,file='mytable1.tab')
Compute
Summary of computational transaction
Raw Input
view raw input (R code)
Raw Output
view raw output of R engine
Computing time
0 seconds
R Server
Big Analytics Cloud Computing Center
Click here to blog (archive) this computation