source('/home/pw/wessanet/cretab')



myrfcuid = 'r0446300'

x <- c(2340,2377.5,2500,2395,2375,2395,2452.5,2462.5,2840,2750,2795,3035,2827.5,2850,3265,3220,3075,3325,3235,3320,3277.5,3050,3432.5,3575,3287.5,3297.5,3255,3135,3177.5,3172.5,3365,3775,3737.5,4060,4197.5,4195,4112.5,3977.5,3982.5,3985,4155,4295,4072.5,4017.5,3965,3847.5,3799.5,4215,4115,4210,4132.5,3955,4057.5,4152.5,4062.5,4092.5,4157.5,4325,4210,5107.5,4862.5,4920,5080,4615,4615,5057.5,4850,4720,5040,4815,5170,6420,5847.5,6177.5,6960,7170,6485,6785,7300,7105,7030,7222.5,7522.5,9500,8697.5,8662.5,8685,8642.5,8492.5,7702.5,7965,8097.5,7947.5,9227.5,7680,7820,7635,7347.5,7177.5,6982.5,7075,7230,7170,7022.5,7087.5,7440,7592.5,8297.5,7852.5,8097.5,8022.5,8062.5,7982.5,9540,8172.5,8007.5)
par3 = 'BFGS'
par2 = '12'
par1 = '12'
par3 <- 'BFGS'
par2 <- '12'
par1 <- '12'
#'GNU S' R Code compiled by R2WASP v. 1.2.327 (Mon, 14 Nov 2016 00:03:48 +0100)
#Author: root
#To cite this work: Wessa P., (2016), Structural Time Series Models (v1.0.6) in Free Statistics Software (v$_version), Office for Research Development and Education, URL http://www.wessa.net/rwasp_structuraltimeseries.wasp/
#Source of accompanying publication: Office for Research, Development, and Education
#
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
postscript(file="/home/pw/wessanet/rcomp/tmp/1cqbf1482326351.ps",horizontal=F,onefile=F,pagecentre=F,paper="special",width=8.3333333333333,height=5.5555555555556) 
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()
postscript(file="/home/pw/wessanet/rcomp/tmp/22nc21482326351.ps",horizontal=F,onefile=F,pagecentre=F,paper="special",width=8.3333333333333,height=5.5555555555556) 
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()
postscript(file="/home/pw/wessanet/rcomp/tmp/38tls1482326351.ps",horizontal=F,onefile=F,pagecentre=F,paper="special",width=8.3333333333333,height=5.5555555555556) 
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()
postscript(file="/home/pw/wessanet/rcomp/tmp/4zx1x1482326351.ps",horizontal=F,onefile=F,pagecentre=F,paper="special",width=8.3333333333333,height=5.5555555555556) 
plot(as.numeric(m$resid),main='Standardized Residuals',ylab='Residuals',xlab='time',type='b')
grid()
dev.off()
postscript(file="/home/pw/wessanet/rcomp/tmp/5qkl01482326351.ps",horizontal=F,onefile=F,pagecentre=F,paper="special",width=8.3333333333333,height=5.5555555555556) 
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()
postscript(file="/home/pw/wessanet/rcomp/tmp/6br681482326351.ps",horizontal=F,onefile=F,pagecentre=F,paper="special",width=8.3333333333333,height=5.5555555555556) 
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()

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="/home/pw/wessanet/rcomp/tmp/74n2y1482326351.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="/home/pw/wessanet/rcomp/tmp/85uqg1482326351.tab") 

try(system("convert /home/pw/wessanet/rcomp/tmp/1cqbf1482326351.ps /home/pw/wessanet/rcomp/tmp/1cqbf1482326351.png",intern=TRUE))
try(system("convert /home/pw/wessanet/rcomp/tmp/22nc21482326351.ps /home/pw/wessanet/rcomp/tmp/22nc21482326351.png",intern=TRUE))
try(system("convert /home/pw/wessanet/rcomp/tmp/38tls1482326351.ps /home/pw/wessanet/rcomp/tmp/38tls1482326351.png",intern=TRUE))
try(system("convert /home/pw/wessanet/rcomp/tmp/4zx1x1482326351.ps /home/pw/wessanet/rcomp/tmp/4zx1x1482326351.png",intern=TRUE))
try(system("convert /home/pw/wessanet/rcomp/tmp/5qkl01482326351.ps /home/pw/wessanet/rcomp/tmp/5qkl01482326351.png",intern=TRUE))
try(system("convert /home/pw/wessanet/rcomp/tmp/6br681482326351.ps /home/pw/wessanet/rcomp/tmp/6br681482326351.png",intern=TRUE))