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(10414.9,12476.8,12384.6,12266.7,12919.9,11497.3,12142,13919.4,12656.8,12034.1,13199.7,10881.3,11301.2,13643.9,12517,13981.1,14275.7,13435,13565.7,16216.3,12970,14079.9,14235,12213.4,12581,14130.4,14210.8,14378.5,13142.8,13714.7,13621.9,15379.8,13306.3,14391.2,14909.9,14025.4,12951.2,14344.3,16093.4,15413.6,14705.7,15972.8,16241.4,16626.4,17136.2,15622.9,18003.9,16136.1,14423.7,16789.4,16782.2,14133.8,12607,12004.5,12175.4,13268,12299.3,11800.6,13873.3,12269.6) > 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 322419.1107 136.6057 156252.8591 0.0000 > m$fitted level slope sea Jan 1 10414.90 0.000000 0.000000 Feb 1 11986.71 72.590769 490.085684 Mar 1 12380.82 88.356635 3.776056 Apr 1 12339.92 84.815177 -73.215338 May 1 12644.25 87.901746 275.646070 Jun 1 11992.21 80.275257 -494.910190 Jul 1 11941.60 78.899920 200.402095 Aug 1 13068.49 90.854193 850.911042 Sep 1 12986.51 88.771723 -329.708844 Oct 1 12382.27 80.161398 -348.173631 Nov 1 12732.66 83.589596 467.040253 Dec 1 11673.17 68.817145 -791.867731 Jan 2 11845.09 66.231335 -543.890758 Feb 2 12708.08 72.335865 935.822435 Mar 2 12741.00 71.330593 -224.004217 Apr 2 13533.68 89.530975 447.423531 May 2 13752.15 91.994881 523.550219 Jun 2 13938.65 93.369742 -503.645451 Jul 2 13924.80 91.996537 -359.097370 Aug 2 14701.08 100.759359 1515.217038 Sep 2 13927.11 89.028697 -957.111373 Oct 2 14125.29 90.522684 -45.393306 Nov 2 13703.07 84.153865 531.925122 Dec 2 13329.35 80.343229 -1115.953563 Jan 3 13385.51 80.268210 -804.506611 Feb 3 13349.48 79.038097 780.916122 Mar 3 14169.69 93.353216 41.111536 Apr 3 14204.85 92.066033 173.650444 May 3 13396.47 73.606952 -253.674872 Jun 3 13804.03 79.575006 -89.333424 Jul 3 14104.30 83.135461 -482.401620 Aug 3 13837.99 77.728237 1541.811373 Sep 3 14057.39 79.894419 -751.093383 Oct 3 14162.38 80.265879 228.823736 Nov 3 14201.52 79.712094 708.376507 Dec 3 14746.52 85.101426 -721.122664 Jan 4 14325.99 79.511882 -1374.787850 Feb 4 14045.08 74.396127 299.223656 Mar 4 15076.17 92.004446 1017.226236 Apr 4 15123.15 91.072923 290.446134 May 4 15158.32 89.918409 -452.617184 Jun 4 15719.94 99.090881 252.861886 Jul 4 16295.45 107.766111 -54.049518 Aug 4 15768.69 96.749164 857.710462 Sep 4 16911.80 114.280870 224.404476 Oct 4 16268.33 102.110496 -645.426673 Nov 4 16858.96 109.519383 1144.936582 Dec 4 16829.40 107.508804 -693.295836 Jan 5 16274.81 97.743480 -1851.107297 Feb 5 16592.44 101.365446 196.959374 Mar 5 16189.03 91.921026 593.174286 Apr 5 14882.59 63.640202 -748.793157 May 5 13877.98 41.666764 -1270.981767 Jun 5 12686.81 16.918397 -682.305830 Jul 5 12214.49 7.484605 -39.086280 Aug 5 12485.16 12.361093 782.835446 Sep 5 12152.46 6.196318 146.843413 Oct 5 12391.65 10.212411 -591.047798 Nov 5 12545.86 12.614520 1327.443377 Dec 5 12653.77 14.180986 -384.169273 > m$resid Jan Feb Mar Apr May Jun 1 0.00000000 1.88245541 0.47712486 -0.22369686 0.38715149 -1.30430891 2 0.19717687 1.38177730 -0.06384193 1.22053363 0.22459373 0.16557734 3 -0.04303126 -0.20221741 1.24768771 -0.09870372 -1.55528192 0.58167140 4 -0.88412628 -0.62380152 1.63101953 -0.07676423 -0.09625666 0.81808710 5 -1.15008811 0.37970502 -0.86471489 -2.39275907 -1.83763748 -2.13245740 Jul Aug Sep Oct Nov Dec 1 -0.23041708 1.84310385 -0.30378726 -1.21762917 0.47463707 -2.00710808 2 -0.18784186 1.19798727 -1.53061689 0.19092572 -0.89583059 -0.80074857 3 0.38506181 -0.60963403 0.24703110 0.04370026 -0.07154174 0.81087785 4 0.82843631 -1.10354942 1.81838504 -1.31530224 0.84748250 -0.24154757 5 -0.84844437 0.45654410 -0.59804372 0.40339842 0.24925418 0.16507282 > 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/1jbbh1260112384.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/2zoaz1260112384.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/3imt81260112384.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/44oxf1260112384.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/5j0ns1260112384.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/69drk1260112384.tab") > system("convert tmp/1jbbh1260112384.ps tmp/1jbbh1260112384.png") > system("convert tmp/2zoaz1260112384.ps tmp/2zoaz1260112384.png") > system("convert tmp/3imt81260112384.ps tmp/3imt81260112384.png") > system("convert tmp/44oxf1260112384.ps tmp/44oxf1260112384.png") > system("convert tmp/5j0ns1260112384.ps tmp/5j0ns1260112384.png") > > > proc.time() user system elapsed 1.491 0.800 2.050