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(21790,13253,37702,30364,32609,30212,29965,28352,25814,22414,20506,28806,22228,13971,36845,35338,35022,34777,26887,23970,22780,17351,21382,24561,17409,11514,31514,27071,29462,26105,22397,23843,21705,18089,20764,25316,17704,15548,28029,29383,36438,32034,22679,24319,18004,17537,20366,22782,19169,13807,29743,25591,29096,26482,22405,27044,17970,18730,19684,19785,18479,10698,31956,29506,34506,27165,26736,23691,18157,17328,18205,20995,17382,9367,31124,26551,30651,25859,25100,25778,20418,18688,20424,24776,19814,12738,31566,30111,30019,31934,25826,26835,20205,17789,20520,22518,15572,11509,25447,24090,27786,26195,20516,22759,19028,16971,20036,22485,18730) > 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 3068393.997 3888.108 3593368.784 0.000 > m$fitted level slope sea Jan 1 21790.00 0.000000 0.00000 Feb 1 16782.01 -165.320321 -3529.01141 Mar 1 26911.88 468.745792 10790.12364 Apr 1 30168.79 608.268779 195.20776 May 1 31906.53 645.157679 702.46840 Jun 1 31465.40 620.787708 -1253.40493 Jul 1 30579.71 592.070795 -614.71274 Aug 1 29248.65 555.108048 -896.65334 Sep 1 27252.65 502.843554 -1438.64778 Oct 1 24446.67 430.834665 -2032.66524 Nov 1 21849.74 361.898664 -1343.73890 Dec 1 24302.83 411.140987 4503.16677 Jan 2 24666.42 411.533990 -2438.42292 Feb 2 24946.95 410.289838 -10975.94637 Mar 2 26511.15 445.289849 10333.84946 Apr 2 31255.03 607.002237 4082.96826 May 2 33507.96 666.369536 1514.04065 Jun 2 34896.75 689.343293 -119.75236 Jul 2 31547.25 574.325083 -4660.24599 Aug 2 27194.77 441.539407 -3224.77189 Sep 2 24038.08 345.798544 -1258.07626 Oct 2 20729.71 249.381958 -3378.70511 Nov 2 21304.03 257.590239 77.96782 Dec 2 20829.82 240.907682 3731.18054 Jan 3 20556.43 230.406179 -3147.42936 Feb 3 22273.52 265.773875 -10759.52293 Mar 3 23625.27 298.203637 7888.72789 Apr 3 24244.85 309.150167 2826.15029 May 3 25935.58 357.808238 3526.42322 Jun 3 25058.44 315.393063 1046.56229 Jul 3 24291.70 279.860067 -1894.70381 Aug 3 24250.24 269.732815 -407.24044 Sep 3 22861.85 219.106985 -1156.85268 Oct 3 22005.07 187.183354 -3916.06816 Nov 3 21092.89 155.545880 -328.88799 Dec 3 21210.69 154.486721 4105.31420 Jan 4 21612.19 161.441119 -3908.19255 Feb 4 24169.56 232.392313 -8621.56268 Mar 4 23158.93 192.808976 4870.06795 Apr 4 24766.20 240.520471 4616.79791 May 4 28396.91 357.969675 8041.08817 Jun 4 30018.52 401.711637 2015.48296 Jul 4 27769.73 311.536742 -5090.73400 Aug 4 25547.74 227.335784 -1228.73877 Sep 4 22002.70 104.864256 -3998.69667 Oct 4 20840.19 64.567914 -3303.18634 Nov 4 20629.06 55.941314 -263.05528 Dec 4 19987.16 34.253101 2794.83873 Jan 5 21763.30 88.805682 -2594.29655 Feb 5 22409.17 106.649322 -8602.16910 Mar 5 24431.40 169.850669 5311.60024 Apr 5 24259.35 158.276473 1331.64810 May 5 22823.38 103.503813 6272.61716 Jun 5 22299.36 81.888156 4182.64430 Jul 5 23485.69 119.671935 -1080.68514 Aug 5 24783.53 159.492983 2260.46905 Sep 5 23466.30 110.255248 -5496.30463 Oct 5 22529.83 75.776093 -3799.82929 Nov 5 21289.05 32.768046 -1605.04809 Dec 5 19730.44 -19.092715 54.56154 Jan 6 20349.57 1.802105 -1870.56956 Feb 6 20393.84 3.207164 -9695.83520 Mar 6 22905.96 87.339934 9050.03806 Apr 6 25333.62 166.764049 4172.38213 May 6 27112.13 221.876460 7393.87150 Jun 6 26022.37 176.935360 1142.62672 Jul 6 26756.22 195.963179 -20.22169 Aug 6 24080.18 98.424716 -389.18034 Sep 6 22965.57 57.523328 -4808.56619 Oct 6 21501.27 6.544530 -4173.27126 Nov 6 20177.95 -37.811286 -1972.95043 Dec 6 20534.43 -24.675116 460.57107 Jan 7 20408.37 -28.060634 -3026.37235 Feb 7 20705.02 -17.162306 -11338.01619 Mar 7 22082.41 29.939633 9041.58868 Apr 7 22747.51 51.515477 3803.49019 May 7 22733.11 49.268065 7917.88627 Jun 7 23403.55 70.473722 2455.44866 Jul 7 23406.35 68.166329 1693.64701 Aug 7 24137.60 90.708915 1640.40208 Sep 7 24275.09 92.293708 -3857.09377 Oct 7 23408.58 59.926118 -4720.58110 Nov 7 22929.18 41.760558 -2505.18425 Dec 7 23539.24 60.890193 1236.76294 Jan 8 23621.91 61.624303 -3807.91160 Feb 8 24274.24 81.579037 -11536.24437 Mar 8 23823.30 63.535667 7742.69722 Apr 8 24911.97 98.366344 5199.03159 May 8 24079.89 66.693673 5939.10744 Jun 8 26377.22 142.674854 5556.78462 Jul 8 25996.82 124.867356 -170.82329 Aug 8 25493.10 103.497590 1341.90319 Sep 8 24529.74 67.294914 -4324.74226 Oct 8 23420.28 27.431233 -5631.27988 Nov 8 23224.13 19.866045 -2704.12985 Dec 8 22407.37 -8.432548 110.62742 Jan 9 21044.73 -54.261656 -5472.73377 Feb 9 21634.89 -32.427766 -10125.89267 Mar 9 20244.19 -78.515090 5202.81156 Apr 9 19343.58 -106.446964 4746.42016 May 9 20731.38 -55.630601 7054.62200 Jun 9 20694.35 -54.997868 5500.65360 Jul 9 20340.74 -65.153891 175.25561 Aug 9 20319.29 -63.668639 2439.71079 Sep 9 21279.34 -28.905490 -2251.33646 Oct 9 21841.59 -8.848294 -4870.59150 Nov 9 21985.32 -3.674754 -1949.31806 Dec 9 21849.57 -8.152198 635.42836 Jan 10 22857.33 26.296733 -4127.32954 > m$resid Jan Feb Mar Apr May Jun 1 0.00000000 -2.26734768 4.27606069 1.44324779 0.63266785 -0.61645841 2 -0.02829651 -0.07553249 0.60450437 2.24916974 0.89328498 0.40051714 3 -0.28817544 0.82592404 0.58752056 0.17207572 0.74731077 -0.67677436 4 0.13629889 1.31595066 -0.67583271 0.76455429 1.83787375 0.68970581 5 0.95508534 0.30469724 1.04281686 -0.18554476 -0.86599712 -0.34206928 6 0.34887588 0.02318314 1.36647044 1.27246837 0.87661112 -0.71468051 7 -0.05533142 0.17710701 0.75969019 0.34567603 -0.03587492 0.33841876 8 0.01187941 0.32203122 -0.29011510 0.55818723 -0.50663875 1.21522178 9 -0.73814341 0.35121126 -0.74000064 -0.44774754 0.81382388 0.01013245 10 0.55361338 Jul Aug Sep Oct Nov Dec 1 -0.85461074 -1.08852210 -1.44067065 -1.86514224 -1.70434521 1.17581739 2 -2.25266536 -2.74896330 -2.00565925 -2.03359957 0.18047058 -0.40665930 3 -0.59706286 -0.17755502 -0.91534072 -0.59288081 -0.60512976 -0.02079864 4 -1.45371716 -1.39140084 -2.06972420 -0.69433380 -0.15097280 -0.38246168 5 0.60383157 0.64470978 -0.80745022 -0.57174438 -0.71899403 -0.86964517 6 0.30397161 -1.56840998 -0.66209477 -0.83006439 -0.72529817 0.21513983 7 -0.03690390 0.36170900 0.02551379 -0.52263079 -0.29396012 0.30984035 8 -0.28511360 -0.34269310 -0.58150421 -0.64124157 -0.12182971 -0.45595879 9 -0.16272309 0.02381576 0.55785887 0.32209349 0.08312802 -0.07196483 10 > 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/1wocb1260550001.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/2dm801260550001.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/36bnw1260550001.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/4e1wl1260550001.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/5cyq11260550001.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/6u8z71260550001.tab") > system("convert tmp/1wocb1260550001.ps tmp/1wocb1260550001.png") > system("convert tmp/2dm801260550001.ps tmp/2dm801260550001.png") > system("convert tmp/36bnw1260550001.ps tmp/36bnw1260550001.png") > system("convert tmp/4e1wl1260550001.ps tmp/4e1wl1260550001.png") > system("convert tmp/5cyq11260550001.ps tmp/5cyq11260550001.png") > > > proc.time() user system elapsed 1.840 0.830 2.104