Send output to:
Browser Blue - Charts White
Browser Black/White
CSV
Data:
4790.92 4795.33 4822.62 4797.52 4822.17 4843.08 4850.79 4827.02 4796.65 4854.96 4870.81 4891.06 4881.38 4921.43 4956.21 4962.81 4949.38 4977.99 4992.73 5009.02 4990.98 5014.96 5022.23 5028.83 4894.36 4918.13 4936.4 4899.87 4862.89 4882.69 4895.46 4883.8 4855.4 4874.33 4880.94 4861.79 4851.44 4840.22 4842.42 4827.02 4749.77 4866.63 4734.37 4726.44 4753.51 4867.29 4793.35 4822.4 4865.09 4987.67 4900.96 4904.71 4889.52 5015.63 4938.81 4924.73 4871.48 4998.24 4891.06 4876.54 4824.15
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
1 seconds
R Server
Big Analytics Cloud Computing Center
Click here to blog (archive) this computation