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(1915,1843,1761,2858,3968,5061,4661,4269,3857,3568,3274,2987,1683,1381,1071,2772,4485,6181,5479,4782,4067,3489,2903,2330,1736,1483,1242,2334,3423,4523,3986,3462,2908,2575,2237,1904,1610,1251,941,2450,3946,5409,4741,4069,3539,3189,2960,2704,1697,1598,1456,2316,3083,4158,3469,2892,2578,2233,1947,2049) > 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.000 362973.080 3396.666 0.000 > m$fitted level slope sea Jan 1 1915.000 0.00000 0.0000000 Feb 1 1843.994 -71.08786 -0.9943178 Mar 1 1760.731 -82.79603 0.2690167 Apr 1 2822.079 1008.51470 35.9214775 May 1 3973.665 1145.05557 -5.6653170 Jun 1 5064.750 1093.53921 -3.7500108 Jul 1 4707.950 -290.80901 -46.9502420 Aug 1 4262.816 -438.11154 6.1842029 Sep 1 3855.577 -408.64356 1.4228863 Oct 1 3565.299 -295.66471 2.7013827 Nov 1 3275.628 -289.94319 -1.6281966 Dec 1 2987.824 -287.90180 -0.8237163 Jan 2 1713.567 -1229.09329 -30.5671224 Feb 2 1310.588 -440.18826 70.4115208 Mar 2 1115.790 -211.95664 -44.7896329 Apr 2 2656.377 1417.39025 115.6227292 May 2 4522.436 1834.19200 -37.4357329 Jun 2 6141.509 1634.28478 39.4906951 Jul 2 5624.430 -364.96892 -145.4298223 Aug 2 4775.853 -814.37350 6.1471592 Sep 2 4056.269 -726.28766 10.7305192 Oct 2 3478.667 -588.11838 10.3331096 Nov 2 2933.570 -548.13916 -30.5700199 Dec 2 2275.398 -650.38812 54.6024138 Jan 3 1789.529 -497.55171 -53.5287803 Feb 3 1358.695 -435.51782 124.3048723 Mar 3 1360.816 -34.35955 -118.8158293 Apr 3 2212.039 779.02828 121.9614138 May 3 3503.317 1248.84228 -80.3169345 Jun 3 4377.375 904.97855 145.6249765 Jul 3 4155.865 -128.62537 -169.8645281 Aug 3 3477.502 -633.02258 -15.5015538 Sep 3 2895.712 -586.01572 12.2880885 Oct 3 2556.156 -359.88405 18.8440511 Nov 3 2254.929 -306.06386 -17.9294447 Dec 3 1868.310 -379.96974 35.6895193 Jan 4 1643.489 -237.65630 -33.4889265 Feb 4 1137.844 -483.58349 113.1555104 Mar 4 1081.907 -94.15839 -140.9069805 Apr 4 2312.197 1114.45403 137.8031812 May 4 4015.432 1651.00214 -69.4320306 Jun 4 5226.227 1249.74236 182.7729479 Jul 4 4952.311 -139.24448 -211.3107199 Aug 4 4111.019 -779.23172 -42.0188268 Sep 4 3523.763 -604.22883 15.2367404 Oct 4 3158.668 -386.23592 30.3321925 Nov 4 2961.765 -213.63958 -1.7651684 Dec 4 2701.527 -256.11333 2.4727919 Jan 5 1728.629 -909.41822 -31.6285932 Feb 5 1424.109 -358.09496 173.8913362 Mar 5 1651.482 173.13618 -195.4821253 Apr 5 2235.674 546.74631 80.3260735 May 5 3160.344 889.87803 -77.3438644 Jun 5 3880.088 735.39116 277.9115024 Jul 5 3665.691 -127.16340 -196.6913993 Aug 5 2968.738 -644.62139 -76.7377498 Sep 5 2556.698 -433.40393 21.3020205 Oct 5 2221.434 -344.27839 11.5662695 Nov 5 1964.376 -265.06936 -17.3764049 Dec 5 1945.147 -41.84546 103.8531656 > m$resid Jan Feb Mar Apr May Jun 1 0.000000000 -0.118782435 -0.019797913 1.809303122 0.226651817 -0.085508159 2 -1.562887978 1.315442868 0.382814720 2.701008189 0.691915328 -0.331810694 3 0.253817081 0.103160041 0.669222723 1.349313696 0.779813126 -0.570766234 4 0.236367686 -0.408410332 0.648005992 2.006024163 0.890443150 -0.666052639 5 -1.085063795 0.915106246 0.882959537 0.620274005 0.569405808 -0.256435030 Jul Aug Sep Oct Nov Dec 1 -2.297778747 -0.244496700 0.048911732 0.187525347 0.009496723 0.003388352 2 -3.318418700 -0.745933776 0.146207240 0.229337083 0.066358629 -0.169715967 3 -1.715601354 -0.837212455 0.078023245 0.375339447 0.089332470 -0.122671608 4 -2.305467174 -1.062268873 0.290474512 0.361830642 0.286481271 -0.070500455 5 -1.431686059 -0.858890513 0.350584462 0.147933213 0.131473533 0.370530377 > 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/15xbu1259871102.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/2xopp1259871102.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/36zgp1259871102.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/4ryr01259871102.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/5d76l1259871102.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/6sq3w1259871102.tab") > > system("convert tmp/15xbu1259871102.ps tmp/15xbu1259871102.png") > system("convert tmp/2xopp1259871102.ps tmp/2xopp1259871102.png") > system("convert tmp/36zgp1259871102.ps tmp/36zgp1259871102.png") > system("convert tmp/4ryr01259871102.ps tmp/4ryr01259871102.png") > system("convert tmp/5d76l1259871102.ps tmp/5d76l1259871102.png") > > > proc.time() user system elapsed 1.300 0.787 1.519