R version 2.13.0 (2011-04-13) Copyright (C) 2011 The R Foundation for Statistical Computing ISBN 3-900051-07-0 Platform: i486-pc-linux-gnu (32-bit) 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(41,39,50,40,43,38,44,35,39,35,29,49,50,59,63,32,39,47,53,60,57,52,70,90,74,62,55,84,94,70,108,139,120,97,126,149,158,124,140,109,114,77,120,133,110,92,97,78,99,107,112,90,98,125,155,190,236,189,174,178,136,161,171,149,184,155,276,224,213,279,268,287) > 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 248.3950263 0.3169035 8.4982038 168.4949790 > m$fitted level slope sea Jan 1 41.00000 0.000000000 0.00000000 Feb 1 39.46818 -0.097929444 -0.12426470 Mar 1 46.54377 0.139008759 0.56204826 Apr 1 42.28031 0.052299543 -0.35396451 May 1 42.67330 0.057832118 0.17583686 Jun 1 39.57156 0.007551655 -0.17063663 Jul 1 42.43876 0.054754542 0.29439199 Aug 1 37.54249 -0.031142262 -0.35220270 Sep 1 38.41981 -0.014570336 0.17890362 Oct 1 36.15323 -0.057637082 -0.15998438 Nov 1 31.36315 -0.152178877 -0.27902344 Dec 1 42.97220 0.092425451 0.85572443 Jan 2 46.51151 -0.066265584 1.93396815 Feb 2 55.71882 0.299433547 0.19120464 Mar 2 60.31502 0.427467212 0.95039115 Apr 2 42.35867 -0.008988013 -2.58617255 May 2 39.52795 -0.071656086 0.67255900 Jun 2 44.51698 0.041188414 0.32961172 Jul 2 49.86579 0.162419590 0.87776471 Aug 2 56.90445 0.323902799 0.17541060 Sep 2 57.03675 0.319278636 0.04452869 Oct 2 53.72863 0.229421798 -0.19151632 Nov 2 64.76143 0.503513431 0.66516734 Dec 2 81.25227 0.902487290 1.98539060 Jan 3 78.98046 0.894980344 -3.62525482 Feb 3 67.94488 0.470231609 -1.56747772 Mar 3 57.88271 0.125994259 1.30522506 Apr 3 75.86995 0.649864970 0.78761732 May 3 88.05637 0.975973464 1.17844993 Jun 3 77.22367 0.642379758 -2.34568071 Jul 3 96.67454 1.179298676 3.56039959 Aug 3 124.70845 1.955658677 3.21184473 Sep 3 122.92025 1.846050103 -1.37656689 Oct 3 107.46268 1.333213932 -3.33210534 Nov 3 119.78377 1.661777083 1.69075897 Dec 3 138.06990 2.144774140 4.09051473 Jan 4 151.69572 2.373433116 1.54626926 Feb 4 134.95792 1.707469621 -3.68439127 Mar 4 138.71154 1.777683332 0.48272049 Apr 4 120.84646 1.139773472 -3.91270594 May 4 114.84038 0.913249464 2.05778257 Jun 4 93.77561 0.219053610 -7.86107446 Jul 4 109.68455 0.716528380 3.95353367 Aug 4 123.25220 1.126293312 4.53868675 Sep 4 115.64435 0.846164694 -2.10508222 Oct 4 103.78168 0.436317881 -6.63314268 Nov 4 99.58748 0.286781077 -0.71232921 Dec 4 83.86844 -0.217539305 0.60332209 Jan 5 88.98441 -0.070006715 7.84236582 Feb 5 103.22050 0.425768674 -1.75891357 Mar 5 107.75445 0.568898366 2.63499431 Apr 5 98.73643 0.244211318 -4.90362600 May 5 95.08028 0.114199945 4.48494652 Jun 5 120.11707 0.941456659 -5.11922946 Jul 5 141.89118 1.633583210 4.74991454 Aug 5 170.26936 2.524310124 9.00135403 Sep 5 214.22897 3.908106135 5.15013802 Oct 5 204.92588 3.465926579 -10.62714894 Nov 5 185.49028 2.700637200 -2.30934728 Dec 5 181.06195 2.466494873 -0.20882555 Jan 6 151.39201 1.463656971 -2.42701385 Feb 6 158.57711 1.661449468 0.18846788 Mar 6 164.17312 1.798695152 5.28663906 Apr 6 157.67266 1.513325093 -5.37419910 May 6 173.25761 1.992312912 5.12978987 Jun 6 167.45692 1.728012195 -9.34851949 Jul 6 235.57674 3.979187215 13.94686110 Aug 6 229.40334 3.634587518 -1.35492719 Sep 6 215.92971 3.053161549 3.89212525 Oct 6 261.08064 4.484799549 1.13481552 Nov 6 268.96432 4.600119557 -2.31860305 Dec 6 280.12092 4.820123150 4.27048225 > m$resid Jan Feb Mar Apr May Jun 1 0.00000000 -0.06389199 0.42890177 -0.27603665 0.02149336 -0.19939485 2 0.26624684 0.49246331 0.25978532 -1.14245177 -0.17589347 0.31538675 3 -0.21691949 -0.67884485 -0.63651373 1.10145767 0.71300634 -0.72964615 4 0.74855994 -1.12123338 0.12362691 -1.20586918 -0.43950016 -1.35141263 5 0.33962795 0.85317549 0.24840184 -0.58717557 -0.23929644 1.52863005 6 -2.01939810 0.34445767 0.23819663 -0.50770181 0.86228710 -0.47738827 Jul Aug Sep Oct Nov Dec 1 0.18032771 -0.31188875 0.05716697 -0.14156051 -0.29716776 0.73778022 2 0.33054030 0.42789605 -0.01191422 -0.22539192 0.67082561 0.99138429 3 1.16161016 1.65776873 -0.23101471 -1.06732218 0.67746299 1.02361044 4 0.96446749 0.78977584 -0.53666059 -0.78075521 -0.28435463 -0.98193785 5 1.27747051 1.63979949 2.54037857 -0.80990376 -1.40323816 -0.43664731 6 4.06610556 -0.62174541 -1.04770783 2.57781715 0.20799794 0.40134418 > 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/wessaorg/rcomp/tmp/1njdi1324143536.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() null device 1 > postscript(file="/var/wessaorg/rcomp/tmp/26vwy1324143536.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() null device 1 > postscript(file="/var/wessaorg/rcomp/tmp/3syh01324143536.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() null device 1 > postscript(file="/var/wessaorg/rcomp/tmp/4r7ow1324143536.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() null device 1 > postscript(file="/var/wessaorg/rcomp/tmp/5bywc1324143536.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() null device 1 > > #Note: the /var/wessaorg/rcomp/createtable file can be downloaded at http://www.wessa.net/cretab > load(file="/var/wessaorg/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/wessaorg/rcomp/tmp/692jq1324143536.tab") > > try(system("convert tmp/1njdi1324143536.ps tmp/1njdi1324143536.png",intern=TRUE)) character(0) > try(system("convert tmp/26vwy1324143536.ps tmp/26vwy1324143536.png",intern=TRUE)) character(0) > try(system("convert tmp/3syh01324143536.ps tmp/3syh01324143536.png",intern=TRUE)) character(0) > try(system("convert tmp/4r7ow1324143536.ps tmp/4r7ow1324143536.png",intern=TRUE)) character(0) > try(system("convert tmp/5bywc1324143536.ps tmp/5bywc1324143536.png",intern=TRUE)) character(0) > > > proc.time() user system elapsed 1.865 0.320 2.190