R version 2.9.0 (2009-04-17) Copyright (C) 2009 The R Foundation for Statistical Computing ISBN 3-900051-07-0 R is free software and comes with ABSOLUTELY NO WARRANTY. You are welcome to redistribute it under certain conditions. Type 'license()' or 'licence()' for distribution details. R is a collaborative project with many contributors. Type 'contributors()' for more information and 'citation()' on how to cite R or R packages in publications. Type 'demo()' for some demos, 'help()' for on-line help, or 'help.start()' for an HTML browser interface to help. Type 'q()' to quit R. > x <- c(0.0314796223103059 + ,-3.00870920563557 + ,-2.07677512619799 + ,-1.25010391965540 + ,0.817975239137125 + ,0.0252076485413113 + ,0.554937772830776 + ,0.230027371950115 + ,2.35672227418686 + ,1.41350455171120 + ,2.73311719024401 + ,1.31551925971717 + ,-2.70076272244080 + ,-0.721411049152714 + ,-0.149388576811997 + ,-0.118199629770334 + ,-0.676562489695275 + ,1.79699928690761 + ,1.79845572032988 + ,0.245100010770855 + ,1.80710848932636 + ,-1.75934771184948 + ,-0.0186697168761931 + ,0.189651523600062 + ,-1.84149562719087 + ,-1.07019530156943 + ,-0.507291477584104 + ,0.866365633831705 + ,-1.76077926699189 + ,-0.580719393339347 + ,-0.435702079860853 + ,-0.994868534845203 + ,1.63136048315789 + ,-1.1949403709466 + ,-1.00525975426991 + ,1.32302234837564 + ,-0.628357549594746 + ,0.632048410440518 + ,-2.16903155809288 + ,2.53779364144266 + ,-0.632933703679292 + ,-1.41749196342200 + ,-0.455343045381255 + ,0.812255211942954 + ,0.627897309219833 + ,0.650904313655623 + ,-1.29800419154382 + ,0.74391671726854 + ,-1.50461634127457 + ,-1.42734677658523 + ,0.263353807408564 + ,-0.430830854870631 + ,0.379576092518008 + ,1.70309353400146 + ,-3.12314448117342 + ,-1.32526207118689 + ,-0.60032490743804 + ,1.23607137604666 + ,0.738007075905376 + ,0.899100896289585) > 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 level slope seas epsilon 0.00000000 0.00000000 0.01397163 1.80429314 > m$fitted level slope sea Jan 1 0.031479622 0.0000000000 0.000000000 Feb 1 -0.337859985 -0.1231132024 -0.153436759 Mar 1 -0.709604204 -0.1852709565 -0.185865736 Apr 1 -0.928920686 -0.1920800617 -0.185290313 May 1 -0.689068229 -0.1200913086 -0.091068338 Jun 1 -0.609751814 -0.0916044909 -0.091244966 Jul 1 -0.425316061 -0.0570994604 -0.044789253 Aug 1 -0.328188352 -0.0399631082 -0.037525313 Sep 1 0.161182279 0.0129702657 0.045591168 Oct 1 0.401884775 0.0336731957 0.033858690 Nov 1 0.833407095 0.0668272894 0.090228037 Dec 1 0.960867064 0.0714913417 0.062453243 Jan 2 0.912510469 0.0629307747 -0.732463639 Feb 2 0.668085692 0.0424404046 0.009096579 Mar 2 0.558253656 0.0329233771 0.030221338 Apr 2 0.476092980 0.0261537269 0.008139327 May 2 0.322020026 0.0161411335 0.013554140 Jun 2 0.536684283 0.0265897190 0.071317353 Jul 2 0.723544337 0.0346032357 0.056216857 Aug 2 0.693826220 0.0315403142 -0.016632352 Sep 2 0.847576942 0.0370953327 0.094889062 Oct 2 0.580412785 0.0238666593 -0.079503984 Nov 2 0.526893172 0.0206422313 0.056468434 Dec 2 0.504624043 0.0189257768 0.034960993 Jan 3 0.389270141 0.0137611738 -0.254724170 Feb 3 0.246067428 0.0079476965 -0.018301178 Mar 3 0.172584824 0.0050394715 0.010430558 Apr 3 0.242548389 0.0072782333 0.049173003 May 3 0.064438908 0.0010986429 -0.114862540 Jun 3 0.001227886 -0.0009758625 0.035187030 Jul 3 -0.041452321 -0.0022791232 0.021239536 Aug 3 -0.119441415 -0.0045733648 -0.093674440 Sep 3 0.010824074 -0.0006075161 0.179733755 Oct 3 -0.075602960 -0.0030595023 -0.171689006 Nov 3 -0.154820058 -0.0051749911 0.017624839 Dec 3 -0.054280831 -0.0023178501 0.127192219 Jan 4 -0.082955349 -0.0030114466 -0.101866706 Feb 4 -0.034454977 -0.0016906307 0.041341888 Mar 4 -0.186727500 -0.0054551781 -0.138095474 Apr 4 -0.011093887 -0.0010383783 0.267429837 May 4 -0.043049620 -0.0017745058 -0.189162105 Jun 4 -0.139152858 -0.0039681972 -0.021717957 Jul 4 -0.167905680 -0.0045314841 0.051616075 Aug 4 -0.103133724 -0.0029914077 -0.057357165 Sep 4 -0.071611530 -0.0022411120 0.202855581 Oct 4 -0.020264583 -0.0011009405 -0.118848634 Nov 4 -0.094724022 -0.0026292426 -0.095723130 Dec 4 -0.058819352 -0.0018428362 0.202286260 Jan 5 -0.119992931 -0.0030294511 -0.208417935 Feb 5 -0.202142046 -0.0045808170 0.041593170 Mar 5 -0.172088230 -0.0039147664 -0.119825199 Apr 5 -0.203879200 -0.0044407325 0.229128886 May 5 -0.169737289 -0.0037262391 -0.095529881 Jun 5 -0.074318741 -0.0019236066 0.085405347 Jul 5 -0.229786100 -0.0046654593 -0.219038677 Aug 5 -0.289519204 -0.0056315584 -0.057348387 Sep 5 -0.320613977 -0.0060705793 0.181570918 Oct 5 -0.247719439 -0.0047321875 0.025823555 Nov 5 -0.198020565 -0.0038250031 -0.088422186 Dec 5 -0.160382719 -0.0031452843 0.256487101 > m$resid Jan Feb Mar Apr May Jun Jul 1 0.0000000 -2.0595575 -0.9950149 -0.1161446 1.3735365 0.6238092 0.8769002 2 -2.1793324 -1.1706906 -0.6148186 -0.4985212 -0.8326268 0.9731528 0.8300655 3 -1.5279429 -1.0452326 -0.5556392 0.4614848 -1.3706142 -0.4935809 -0.3317095 4 -0.3456121 0.4973277 -1.4675970 1.8136814 -0.3182283 -0.9969479 -0.2687406 5 -0.9202491 -1.0033433 0.4400156 -0.3612183 0.5104105 1.3384862 -2.1143771 Aug Sep Oct Nov Dec 1 0.5066637 1.8165915 0.8207363 1.5093011 0.2422822 2 -0.3507038 0.6993237 -1.8223560 -0.4840136 -0.2805034 3 -0.6230853 1.1466103 -0.7530760 -0.6889079 0.9903798 4 0.7703432 0.3929869 0.6246258 -0.8750168 0.4737984 5 -0.7731270 -0.3643180 1.1509268 0.8082994 0.6329714 > 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/1qofq1259788745.ps",horizontal=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() null device 1 > postscript(file="/var/www/html/rcomp/tmp/2jhmi1259788745.ps",horizontal=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() null device 1 > postscript(file="/var/www/html/rcomp/tmp/36hpb1259788745.ps",horizontal=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() null device 1 > postscript(file="/var/www/html/rcomp/tmp/4aj9r1259788745.ps",horizontal=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() null device 1 > postscript(file="/var/www/html/rcomp/tmp/591151259788745.ps",horizontal=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() null device 1 > > #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/6gmu91259788745.tab") > system("convert tmp/1qofq1259788745.ps tmp/1qofq1259788745.png") > system("convert tmp/2jhmi1259788745.ps tmp/2jhmi1259788745.png") > system("convert tmp/36hpb1259788745.ps tmp/36hpb1259788745.png") > system("convert tmp/4aj9r1259788745.ps tmp/4aj9r1259788745.png") > system("convert tmp/591151259788745.ps tmp/591151259788745.png") > > > proc.time() user system elapsed 1.427 0.809 1.757