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(3111,3995,5245,5588,10681,10516,7496,9935,10249,6271,3616,3724,2886,3318,4166,6401,9209,9820,7470,8207,9564,5309,3385,3706,2733,3045,3449,5542,10072,9418,7516,7840,10081,4956,3641,3970,2931,3170,3889,4850,8037,12370,6712,7297,10613,5184,3506,3810,2692,3073,3713,4555,7807,10869,9682,7704,9826,5456,3677,3431,2765,3483,3445,6081,8767,9407,6551,12480,9530,5960,3252,3717,2642,2989,3607,5366,8898,9435,7328,8594,11349,5797,3621,3851) > 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.0 1107282.7 196419.6 557276.7 > m$fitted level slope sea Jan 1 3111.000 0.000000 0.00000 Feb 1 3612.883 488.485602 139.64198 Mar 1 5006.457 1081.271150 89.10848 Apr 1 5743.849 862.125117 -99.94211 May 1 9642.940 2812.721117 559.03180 Jun 1 11121.047 1955.496568 -396.73404 Jul 1 8825.724 -777.585403 -665.40490 Aug 1 9281.956 16.298586 460.07592 Sep 1 10056.153 503.999510 74.31628 Oct 1 7365.667 -1551.396598 -595.13220 Nov 1 4025.625 -2702.156114 -129.94902 Dec 1 3044.826 -1594.716093 410.02649 Jan 2 2486.992 -943.350248 240.71244 Feb 2 2905.409 -135.699925 247.00709 Mar 2 3737.895 444.853436 296.45689 Apr 2 6377.463 1806.628252 -285.13136 May 2 8490.685 1996.402142 675.79585 Jun 2 9763.738 1549.842976 157.02198 Jul 2 9186.053 234.051729 -1418.05247 Aug 2 8276.983 -473.931037 90.39210 Sep 2 8563.549 -2.751793 893.72902 Oct 2 6381.095 -1353.072875 -766.23080 Nov 2 3992.717 -1994.031379 -462.53688 Dec 2 2985.204 -1384.491822 582.70310 Jan 3 2482.907 -839.434634 125.32600 Feb 3 2683.535 -213.210725 226.30355 Mar 3 3405.599 346.098658 -80.46201 Apr 3 5401.355 1352.492519 -84.14865 May 3 8652.616 2513.150791 1162.84893 Jun 3 9541.190 1523.116415 95.98350 Jul 3 9353.042 479.622380 -1605.19276 Aug 3 8410.247 -389.062977 -377.10425 Sep 3 8419.588 -145.597512 1607.28687 Oct 3 6223.908 -1398.150910 -989.43805 Nov 3 4325.725 -1703.327133 -616.88095 Dec 3 3238.198 -1327.964543 648.26389 Jan 4 2703.314 -843.190162 119.21575 Feb 4 2829.666 -258.374033 212.58730 Mar 4 3983.961 587.470541 -281.21137 Apr 4 5264.871 1007.321039 -508.04438 May 4 6489.296 1139.221108 1518.62943 Jun 4 10742.692 3026.787240 1211.70976 Jul 4 9759.169 596.752782 -2510.72087 Aug 4 8350.230 -620.077919 -784.32444 Sep 4 8200.494 -334.523121 2349.40124 Oct 4 6521.232 -1150.726491 -1156.85866 Nov 4 4464.298 -1700.208855 -836.85079 Dec 4 3100.715 -1496.136581 664.11871 Jan 5 2453.835 -980.292615 123.83260 Feb 5 2725.902 -225.489507 181.68539 Mar 5 3750.943 524.035897 -202.73161 Apr 5 5079.658 1009.398696 -631.97912 May 5 6832.955 1459.667489 874.86824 Jun 5 8613.179 1653.409498 2213.23463 Jul 5 11385.657 2329.245344 -1852.43951 Aug 5 10059.319 119.888134 -1868.38794 Sep 5 7831.793 -1299.717723 2307.13586 Oct 5 6306.450 -1436.133641 -820.37855 Nov 5 4602.276 -1598.091259 -889.56209 Dec 5 2917.083 -1650.741617 525.54056 Jan 6 2439.030 -941.204749 169.38971 Feb 6 3036.257 -15.067894 243.58513 Mar 6 3670.202 374.101697 -310.72136 Apr 6 6194.043 1668.037853 -398.53477 May 6 8163.635 1850.101974 563.27977 Jun 6 8214.342 764.557270 1431.15647 Jul 6 7852.620 85.745352 -1152.38355 Aug 6 12141.255 2619.538416 -218.87035 Sep 6 9351.014 -643.280447 897.27901 Oct 6 7006.639 -1669.134089 -820.76999 Nov 6 4349.177 -2265.008163 -965.92127 Dec 6 3029.042 -1694.974079 562.30616 Jan 7 2429.107 -1034.125721 67.37485 Feb 7 2565.099 -330.488692 269.50125 Mar 7 3909.192 673.674180 -522.80671 Apr 7 5516.295 1234.746185 -273.92446 May 7 7618.760 1757.663350 1164.14633 Jun 7 8432.837 1189.201350 1127.04302 Jul 7 9519.740 1127.632916 -2178.21376 Aug 7 8327.309 -268.547959 573.67204 Sep 7 9351.004 509.289307 1826.89503 Oct 7 7312.851 -1024.073334 -1178.46392 Nov 7 5075.713 -1754.312161 -1293.96450 Dec 7 3410.208 -1700.815465 429.01148 > m$resid Jan Feb Mar Apr May Jun 1 0.00000000 0.62018954 0.57252707 -0.21528110 1.86347005 -0.81436023 2 0.67855720 0.72606543 0.56071053 1.29363634 0.17886942 -0.42340515 3 0.53293447 0.58336351 0.53485770 0.96037214 1.09755123 -0.93819434 4 0.46538785 0.55084947 0.80579458 0.40075197 0.12504257 1.78898382 5 0.49209311 0.71390086 0.71302162 0.46293035 0.42752091 0.18371013 6 0.67520886 0.87752637 0.37001253 1.23314968 0.17301477 -1.02988109 7 0.62821948 0.66730529 0.95451709 0.53438216 0.49714818 -0.53956334 Jul Aug Sep Oct Nov Dec 1 -2.59759688 0.75451060 0.46347630 -1.95328657 -1.09359249 1.05242471 2 -1.25048211 -0.67282797 0.44775750 -1.28325385 -0.60911729 0.58013293 3 -0.99141937 -0.82557973 0.23136489 -1.19033085 -0.29006500 0.35771117 4 -2.30796268 -1.15641550 0.27136911 -0.77566669 -0.52243825 0.19446728 5 0.64172293 -2.09946330 -1.34911384 -0.12965174 -0.15402532 -0.05014509 6 -0.64446803 2.40746254 -3.10085236 -0.97511085 -0.56676549 0.54260060 7 -0.05845182 -1.32641208 0.73923376 -1.45770337 -0.69459643 0.05089865 > 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/1jme91292749221.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/www/html/rcomp/tmp/2cvvc1292749221.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/www/html/rcomp/tmp/3cvvc1292749221.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/www/html/rcomp/tmp/45muf1292749221.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/www/html/rcomp/tmp/55muf1292749221.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/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/61wa51292749221.tab") > > try(system("convert tmp/1jme91292749221.ps tmp/1jme91292749221.png",intern=TRUE)) character(0) > try(system("convert tmp/2cvvc1292749221.ps tmp/2cvvc1292749221.png",intern=TRUE)) character(0) > try(system("convert tmp/3cvvc1292749221.ps tmp/3cvvc1292749221.png",intern=TRUE)) character(0) > try(system("convert tmp/45muf1292749221.ps tmp/45muf1292749221.png",intern=TRUE)) character(0) > try(system("convert tmp/55muf1292749221.ps tmp/55muf1292749221.png",intern=TRUE)) character(0) > > > proc.time() user system elapsed 1.542 0.857 3.771