x <- c(19876,45335,48674,156392,100837,101605,532850,294189,80763,105995,25045,90474,48481,50730,68694,207716,99132,104012,422632,364974,82687,66834,28408,97073,40284,24421,116346,72120,108751,91738,402216,390070,106045,110070,70668,167841,28607,95371,30605,131063,81214,85451,455196,454570,63114,74287,42350,113375)
par1 = '12'
#'GNU S' R Code compiled by R2WASP v. 1.0.44 ()
#Author: Prof. Dr. P. Wessa
#To cite this work: AUTHOR(S), (YEAR), YOUR SOFTWARE TITLE (vNUMBER) in Free Statistics Software (v$_version), Office for Research Development and Education, URL http://www.wessa.net/rwasp_YOURPAGE.wasp/
#Source of accompanying publication: Office for Research, Development, and Education
#Technical description: Write here your technical program description (don't use hard returns!)
par1 <- as.numeric(par1)
nx <- length(x)
x <- ts(x,frequency=par1)
m <- StructTS(x,type='BSM')
m$coef
m$fitted
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
mylagmax <- nx/2
postscript(file="/var/www/html/rcomp/tmp/1w1xh1292764750.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="/var/www/html/rcomp/tmp/2w1xh1292764750.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="/var/www/html/rcomp/tmp/3psw11292764750.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="/var/www/html/rcomp/tmp/4psw11292764750.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="/var/www/html/rcomp/tmp/502wm1292764750.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()

#Note: the /var/www/html/rcomp/createtable file can be downloaded at http://www.wessa.net/cretab
load(file="/var/www/html/rcomp/createtable")

a<-table.start()
a<-table.row.start(a)
a<-table.element(a,'Structural Time Series Model',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="/var/www/html/rcomp/tmp/6etcv1292764750.tab") 

try(system("convert tmp/1w1xh1292764750.ps tmp/1w1xh1292764750.png",intern=TRUE))
try(system("convert tmp/2w1xh1292764750.ps tmp/2w1xh1292764750.png",intern=TRUE))
try(system("convert tmp/3psw11292764750.ps tmp/3psw11292764750.png",intern=TRUE))
try(system("convert tmp/4psw11292764750.ps tmp/4psw11292764750.png",intern=TRUE))
try(system("convert tmp/502wm1292764750.ps tmp/502wm1292764750.png",intern=TRUE))

