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



myrfcuid = 'r0657140'

x <- c(3370,2740,2920,2730,2920,4260,3210,6280,6930,4190,4740,4800,3850,2970,3450,2150,1900,2910,2200,3570,5120,4980,3990,3350,3850,3520,1790,1660,1540,2210,2060,2580,4590,4120,4700,3470,2850,1950,3190,2130,1350,1510,1470,3210,2150,3210,2540,2330,3230,1810,2060,1360,990,2310,1020,1350,2570,3170,2170,1370,1330,770,2160,1520,1070,810,1030,2640,2060,4180,2160,1500,1930,2300,1680,1070,700,1390,1250,1120,2170,1510,1630,960,1700,1390,880,860,800,1410,1200,1240,1920,1440,1940,1380,1660,1450,820,1020,820,1320,780,1340,1630,2210,2010,2100)
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/1mh861482238737.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/2lh6m1482238737.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/3p4iq1482238737.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/4lwg11482238737.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/5bz771482238737.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/6xxb11482238737.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/7sr9p1482238737.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/8ya111482238737.tab") 

try(system("convert /home/pw/wessanet/rcomp/tmp/1mh861482238737.ps /home/pw/wessanet/rcomp/tmp/1mh861482238737.png",intern=TRUE))
try(system("convert /home/pw/wessanet/rcomp/tmp/2lh6m1482238737.ps /home/pw/wessanet/rcomp/tmp/2lh6m1482238737.png",intern=TRUE))
try(system("convert /home/pw/wessanet/rcomp/tmp/3p4iq1482238737.ps /home/pw/wessanet/rcomp/tmp/3p4iq1482238737.png",intern=TRUE))
try(system("convert /home/pw/wessanet/rcomp/tmp/4lwg11482238737.ps /home/pw/wessanet/rcomp/tmp/4lwg11482238737.png",intern=TRUE))
try(system("convert /home/pw/wessanet/rcomp/tmp/5bz771482238737.ps /home/pw/wessanet/rcomp/tmp/5bz771482238737.png",intern=TRUE))
try(system("convert /home/pw/wessanet/rcomp/tmp/6xxb11482238737.ps /home/pw/wessanet/rcomp/tmp/6xxb11482238737.png",intern=TRUE))