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(49915,47469,45652,43492,41087,42931,67256,72316,65624,59450,52851,51214,44092,43752,40320,40551,38329,39530,59648,61031,55560,43877,38510,36085,35994,32617,30001,27894,26083,28817,48742,49915,40264,34276,30426,30793,29855,28081,26820,25782,22654,27373,43675,45096,38145,34017,31537,33814,36531,36935,36497,35110,33137,37407,53963,56602,49694,43957,41723,45599,42503,42153,39098,37449,34748,36548,53639,55289,47774,42156,38019) > par20 = '' > par19 = '' > par18 = '' > par17 = '' > par16 = '' > par15 = '' > par14 = '' > par13 = '' > par12 = '' > par11 = '' > par10 = '' > par9 = '' > par8 = '' > par7 = '' > par6 = '' > par5 = '' > par4 = '' > par3 = '' > par2 = '' > par1 = '12' > ylab = '' > xlab = '' > main = '' > #'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 12083351 2452157 0 > m$fitted level slope sea Jan 1 49915.00 0.00000 0.00000 Feb 1 48054.91 -1908.51229 -585.91003 Mar 1 45445.84 -2397.35531 206.16437 Apr 1 43332.34 -2187.19967 159.65645 May 1 41106.68 -2215.08495 -19.67513 Jun 1 42052.44 64.93651 878.56023 Jul 1 62050.99 14488.87330 5205.01487 Aug 1 74531.93 13035.11297 -2215.93487 Sep 1 70219.95 477.10442 -4595.94963 Oct 1 60738.65 -6730.69754 -1288.64567 Nov 1 52451.88 -7856.89421 399.12393 Dec 1 49691.56 -4168.21377 1522.44458 Jan 2 44728.48 -4742.06001 -636.47501 Feb 2 43191.99 -2419.20896 560.01483 Mar 2 40188.52 -2827.23953 131.47990 Apr 2 38604.22 -1960.67443 1946.77744 May 2 38299.04 -799.62244 29.95959 Jun 2 43672.16 3495.66131 -4142.16253 Jul 2 53673.29 8021.82547 5974.70771 Aug 2 59842.37 6730.77302 1188.62955 Sep 2 59281.74 1648.19522 -3721.74467 Oct 2 47941.10 -7405.47819 -4064.09552 Nov 2 39063.84 -8430.99907 -553.83623 Dec 2 33352.91 -6537.60563 2732.09265 Jan 3 35089.75 -777.40712 904.25438 Feb 3 32448.58 -2076.48225 168.41774 Mar 3 29617.93 -2596.69678 383.07458 Apr 3 26136.28 -3205.52354 1757.72231 May 3 26954.11 -423.97716 -871.11042 Jun 3 33937.80 4682.14216 -5120.79931 Jul 3 42505.88 7355.45030 6236.11833 Aug 3 47606.19 5802.45734 2308.80926 Sep 3 42698.50 -1578.18059 -2434.49526 Oct 3 37913.02 -3788.44063 -3637.01606 Nov 3 32318.95 -5032.14095 -1892.94734 Dec 3 29861.03 -3260.16340 931.96780 Jan 4 27944.69 -2334.47521 1910.31294 Feb 4 26763.08 -1540.58725 1317.91637 Mar 4 25427.34 -1400.23322 1392.65571 Apr 4 24835.06 -847.44908 946.94351 May 4 25819.37 409.44069 -3165.36573 Jun 4 32765.20 4891.25097 -5392.19621 Jul 4 37433.73 4738.80996 6241.27436 Aug 4 40069.84 3299.28662 5026.16132 Sep 4 40244.81 1158.92368 -2099.80765 Oct 4 37902.81 -1239.66627 -3885.80991 Nov 4 34644.36 -2622.29164 -3107.36266 Dec 4 32934.16 -1997.71739 879.84267 Jan 5 33968.27 79.66295 2562.73181 Feb 5 35099.43 799.63445 1835.57310 Mar 5 35141.91 282.73763 1355.08543 Apr 5 34899.79 -75.28951 210.20994 May 5 37795.49 1954.96631 -4658.49329 Jun 5 42466.89 3811.37294 -5059.88701 Jul 5 47190.50 4434.01297 6772.49833 Aug 5 50896.44 3937.19309 5705.56275 Sep 5 51468.28 1639.83201 -1774.27708 Oct 5 48267.75 -1664.88037 -4310.74707 Nov 5 45313.13 -2545.29819 -3590.13183 Dec 5 45066.03 -976.29330 532.97287 Jan 6 41079.21 -3032.42102 1423.78510 Feb 6 39241.97 -2216.74050 2911.02583 Mar 6 37323.88 -2013.31672 1774.12085 Apr 6 37534.09 -499.51286 -85.08642 May 6 39775.05 1368.48647 -5027.05375 Jun 6 42079.27 2006.50762 -5531.26990 Jul 6 46329.90 3535.34721 7309.09839 Aug 6 48919.24 2891.17195 6369.75629 Sep 6 48789.95 834.15109 -1015.94635 Oct 6 46880.64 -1034.44293 -4724.64150 Nov 6 42992.60 -2978.10310 -4973.60127 > m$resid Jan Feb Mar Apr May Jun 1 0.000000000 -0.625488586 -0.149024277 0.062131387 -0.007956513 0.656803020 2 -0.165575794 0.678694016 -0.116793935 0.253193319 0.333835640 1.232966714 3 1.662039712 -0.374255852 -0.149313969 -0.176114056 0.801730507 1.466135027 4 0.266779262 0.228354022 0.040331477 0.159424303 0.362209382 1.288049911 5 0.598185647 0.207037952 -0.148605247 -0.103140945 0.584835371 0.533830024 6 -0.591783984 0.234553622 0.058497052 0.435879677 0.537931762 0.183526345 Jul Aug Sep Oct Nov Dec 1 4.151212701 -0.418188710 -3.612739860 -2.073529237 -0.323981539 1.061151691 2 1.303063376 -0.371467633 -1.462051070 -2.604713678 -0.295022140 0.544858643 3 0.768965230 -0.446898220 -2.123177000 -0.635865599 -0.357789863 0.510085539 4 -0.043829764 -0.414172084 -0.615752764 -0.690023241 -0.397788230 0.179793888 5 0.179012463 -0.142913943 -0.660928068 -0.950733352 -0.253325187 0.451620286 6 0.439587940 -0.185278399 -0.591773604 -0.537612855 -0.559299923 > 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/18dsv1293564363.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/28dsv1293564363.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/3j5rg1293564363.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/4j5rg1293564363.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/5ce8j1293564363.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/6xw7p1293564363.tab") > > try(system("convert tmp/18dsv1293564363.ps tmp/18dsv1293564363.png",intern=TRUE)) character(0) > try(system("convert tmp/28dsv1293564363.ps tmp/28dsv1293564363.png",intern=TRUE)) character(0) > try(system("convert tmp/3j5rg1293564363.ps tmp/3j5rg1293564363.png",intern=TRUE)) character(0) > try(system("convert tmp/4j5rg1293564363.ps tmp/4j5rg1293564363.png",intern=TRUE)) character(0) > try(system("convert tmp/5ce8j1293564363.ps tmp/5ce8j1293564363.png",intern=TRUE)) character(0) > > > proc.time() user system elapsed 1.330 0.856 5.458