R version 2.12.0 (2010-10-15) Copyright (C) 2010 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(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,14538) > 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 1525416 2125017 1797465 0 > m$fitted level slope sea Jan 1 31514.00 0.00000 0.000000 Feb 1 28488.22 -1092.63512 -1417.221227 Mar 1 28256.63 -703.14339 1205.365028 Apr 1 26681.19 -1151.42398 -576.186345 May 1 23266.15 -2311.89232 -869.145656 Jun 1 22753.24 -1398.55116 1089.757499 Jul 1 21809.55 -1167.00062 -104.551961 Aug 1 18952.92 -2029.74343 -863.915487 Sep 1 19423.46 -752.00922 1340.541755 Oct 1 23548.10 1739.87741 1767.902472 Nov 1 20588.72 -660.82119 -2884.716631 Dec 1 16432.92 -2446.08986 -884.915899 Jan 2 21943.66 1543.87382 6085.341909 Feb 2 29386.94 4496.42367 -3.943106 Mar 2 35122.86 5090.31100 1315.136394 Apr 2 34834.21 2469.73347 -2800.212377 May 2 28115.16 -2084.13305 -5436.158912 Jun 2 23375.79 -3394.70452 943.205001 Jul 2 18142.40 -4299.60362 -138.400328 Aug 2 17756.71 -2369.74896 -219.707218 Sep 2 19234.97 -468.60728 1131.032215 Oct 2 18961.24 -372.29044 3820.760285 Nov 2 20082.62 364.86635 -913.620800 Dec 2 19011.61 -342.28322 -5204.612448 Jan 3 23961.36 2270.02752 5781.639197 Feb 3 26781.71 2540.76813 -1190.712499 Mar 3 26651.44 1247.05485 2444.564979 Apr 3 26083.70 367.49364 398.302101 May 3 26319.78 303.19325 -3914.776979 Jun 3 25817.57 -91.22784 1226.425758 Jul 3 21602.58 -2104.92348 -3632.578062 Aug 3 19759.92 -1976.89123 -1029.922527 Sep 3 18318.03 -1715.42209 1365.970287 Oct 3 16370.96 -1828.66142 3414.037099 Nov 3 17692.56 -289.87940 786.437912 Dec 3 18138.72 69.69605 -7440.715016 Jan 4 23764.08 2788.24887 8191.915872 Feb 4 29092.15 4027.60057 413.850679 Mar 4 32248.88 3605.67405 2257.119083 Apr 4 29681.18 617.43351 -2516.177024 May 4 29778.33 364.35027 -3042.325864 Jun 4 23811.59 -2721.31135 -120.592199 Jul 4 21100.11 -2716.52817 -2943.111590 Aug 4 18284.83 -2764.51201 -956.827837 Sep 4 16196.06 -2436.20072 2008.936912 Oct 4 17041.48 -841.63714 3953.518295 Nov 4 17150.35 -379.71352 231.653173 Dec 4 18676.40 547.31698 -9309.396533 Jan 5 22814.11 2295.14738 8309.888325 Feb 5 25761.01 2611.71108 789.989703 Mar 5 26900.78 1899.69548 3750.218797 Apr 5 27990.59 1508.24232 -2131.587982 May 5 26773.81 187.37900 -1673.807453 Jun 5 25929.12 -313.81095 -151.119048 Jul 5 23857.14 -1166.84705 -3439.135746 Aug 5 20657.79 -2151.35322 -1969.792107 Sep 5 19060.19 -1883.35131 1363.809902 Oct 5 19961.14 -535.81837 4814.858198 Nov 5 20438.03 -45.36594 -624.026705 Dec 5 22675.84 1061.62281 -9937.835454 Jan 6 23795.62 1089.83281 7770.376983 Feb 6 27632.20 2419.55305 2478.795370 Mar 6 27424.71 1151.24044 2594.289180 Apr 6 31443.63 2534.91515 490.368663 May 6 29391.15 317.33513 -3565.152544 Jun 6 26809.44 -1086.28065 25.556929 Jul 6 23458.15 -2182.41666 -3253.147857 Aug 6 20273.89 -2666.48518 -2484.889893 Sep 6 19320.77 -1839.58316 1199.226657 Oct 6 18071.61 -1554.66205 4446.388081 Nov 6 16923.68 -1358.21542 -1351.683535 Dec 6 19836.25 706.78714 -8327.246046 Jan 7 19628.27 264.42911 5818.731214 Feb 7 20637.07 623.83743 3452.931891 Mar 7 24742.26 2301.42089 3043.739295 Apr 7 25022.48 1327.87850 1172.516457 May 7 23679.83 39.84567 -3163.832460 Jun 7 22113.10 -736.03069 645.901829 Jul 7 21634.76 -611.61216 -2606.758106 Aug 7 20354.39 -934.05159 -3383.394968 Sep 7 19036.55 -1118.87350 999.447907 Oct 7 18251.63 -958.07479 4233.365343 Nov 7 20300.45 490.97072 -1570.447471 Dec 7 22239.33 1189.43282 -7701.331868 > m$resid Jan Feb Mar Apr May Jun 1 0.000000000 -1.465983005 0.311490301 -0.309348729 -0.773902918 0.621016281 2 2.886379570 1.997365804 0.403221519 -1.832777179 -3.122368285 -0.891660310 3 1.810868667 0.185062379 -0.882823073 -0.608103735 -0.044302416 -0.269655335 4 1.870667684 0.848002732 -0.288666585 -2.057354680 -0.174302631 -2.115445645 5 1.200075751 0.216709054 -0.487645386 -0.269095944 -0.908922851 -0.343984938 6 0.019353756 0.910600280 -0.869069542 0.950530719 -1.524966521 -0.963697358 7 -0.303393891 0.246192741 1.149822692 -0.668539731 -0.885312009 -0.532730012 Jul Aug Sep Oct Nov Dec 1 0.158797926 -0.591919735 0.876493082 1.709405567 -1.646861834 -1.224678577 2 -0.619028940 1.323840044 1.304212227 0.066072296 0.505788839 -0.486851166 3 -1.376823448 0.087748466 0.179359382 -0.077695110 1.056707119 0.247479032 4 0.003273047 -0.032867975 0.225156364 1.094440704 0.317369819 0.637433920 5 -0.584207483 -0.674280898 0.183761235 0.925086774 0.337018950 0.760683058 6 -0.751122983 -0.331567149 0.566935581 0.195608480 0.134982435 1.418427184 7 0.085287105 -0.220897776 -0.126715898 0.110388169 0.995514662 0.479639833 > 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/rcomp/tmp/1xlu01291912029.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/rcomp/tmp/28ccl1291912029.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/rcomp/tmp/313to1291912029.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/rcomp/tmp/413to1291912029.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/rcomp/tmp/513to1291912029.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/rcomp/createtable file can be downloaded at http://www.wessa.net/cretab > load(file="/var/www/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/rcomp/tmp/6xd9x1291912029.tab") > > try(system("convert tmp/1xlu01291912029.ps tmp/1xlu01291912029.png",intern=TRUE)) character(0) > try(system("convert tmp/28ccl1291912029.ps tmp/28ccl1291912029.png",intern=TRUE)) character(0) > try(system("convert tmp/313to1291912029.ps tmp/313to1291912029.png",intern=TRUE)) character(0) > try(system("convert tmp/413to1291912029.ps tmp/413to1291912029.png",intern=TRUE)) character(0) > try(system("convert tmp/513to1291912029.ps tmp/513to1291912029.png",intern=TRUE)) character(0) > > > proc.time() user system elapsed 1.580 0.930 2.475