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(4716.99,4926.65,4920.10,5170.09,5246.24,5283.61,4979.05,4825.20,4695.12,4711.54,4727.22,4384.96,4378.75,4472.93,4564.07,4310.54,4171.38,4049.38,3591.37,3720.46,4107.23,4101.71,4162.34,4136.22,4125.88,4031.48,3761.36,3408.56,3228.47,3090.45,2741.14,2980.44,3104.33,3181.57,2863.86,2898.01,3112.33,3254.33,3513.47,3587.61,3727.45,3793.34,3817.58,3845.13,3931.86,4197.52,4307.13,4229.43,4362.28,4217.34,4361.28,4327.74,4417.65,4557.68,4650.35,4967.18,5123.42,5290.85,5535.66,5514.06,5493.88,5694.83,5850.41,6116.64,6175.00,6513.58,6383.78,6673.66,6936.61,7300.68,7392.93,7497.31,7584.71,7160.79,7196.19,7245.63,7347.51,7425.75,7778.51,7822.33,8181.22,8371.47,8347.71,8672.11,8802.79,9138.46,9123.29,9023.21,8850.41,8864.58,9163.74,8516.66,8553.44,7555.20,7851.22,7442.00,7992.53,8264.04,7517.39,7200.40,7193.69,6193.58,5104.21,4800.46,4461.61,4398.59,4243.63,4293.82) > 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 62908.665 1366.889 0.000 0.000 > m$fitted level slope sea Jan 1 4716.990 0.000000 0.0000000 Feb 1 4915.824 12.291777 10.8261065 Mar 1 4909.327 11.478358 10.7731667 Apr 1 5158.688 26.014198 11.4024515 May 1 5234.715 29.842534 11.5246282 Jun 1 5272.069 30.514930 11.5413335 Jul 1 4968.178 -2.977031 10.8720493 Aug 1 4814.597 -19.343967 10.6033847 Sep 1 4684.691 -32.101695 10.4289115 Oct 1 4701.044 -26.257573 10.4961533 Nov 1 4716.673 -21.037399 10.5470366 Dec 1 4374.753 -62.010691 10.2070069 Jan 2 4436.733 -47.704126 -57.9834387 Feb 2 4465.896 -37.412346 7.0341641 Mar 2 4557.042 -20.331740 7.0284906 Apr 2 4303.503 -51.556277 7.0373987 May 2 4164.340 -63.353649 7.0402945 Jun 2 4042.338 -71.285219 7.0419708 Jul 2 3584.318 -123.754103 7.0515255 Aug 2 3713.414 -89.368108 7.0461281 Sep 2 4100.193 -24.500480 7.0373489 Oct 2 4094.673 -21.911236 7.0370467 Nov 2 4155.304 -10.640252 7.0359119 Dec 2 4129.184 -12.755549 7.0360956 Jan 3 4188.239 -3.116789 -62.3589375 Feb 3 4028.101 -24.381011 3.3792952 Mar 3 3757.921 -58.014836 3.4388135 Apr 3 3405.060 -98.366960 3.5004368 May 3 3224.955 -109.554868 3.5151818 Jun 3 3086.930 -113.452045 3.5196145 Jul 3 2737.589 -145.745287 3.5513145 Aug 3 2976.933 -93.023398 3.5066493 Sep 3 3100.845 -63.321859 3.4849328 Oct 3 3178.097 -44.074547 3.4727872 Nov 3 2860.367 -81.544436 3.4931936 Dec 3 2894.524 -65.701777 3.4857471 Jan 4 3111.681 -27.292940 0.6492892 Feb 4 3252.582 -4.491300 1.7480195 Mar 4 3511.776 31.621992 1.6935441 Apr 4 3585.924 37.445809 1.6859631 May 4 3725.780 51.470195 1.6702067 Jun 4 3791.672 53.445109 1.6682917 Jul 4 3815.908 49.445356 1.6716391 Aug 4 3843.456 46.446777 1.6738050 Sep 4 3930.190 51.963480 1.6703659 Oct 4 4195.865 81.228366 1.6546204 Nov 4 4305.477 85.115078 1.6528156 Dec 4 4227.768 62.818604 1.6617513 Jan 5 4372.183 73.927908 -9.9026590 Feb 5 4218.379 43.000632 -1.0389371 Mar 5 4362.336 56.827256 -1.0558791 Apr 5 4328.783 44.449695 -1.0427895 May 5 4418.698 50.676035 -1.0484725 Jun 5 4558.738 62.913649 -1.0581128 Jul 5 4651.411 66.988855 -1.0608835 Aug 5 4968.261 101.204454 -1.0809612 Sep 5 5124.505 108.741414 -1.0847782 Oct 5 5291.938 116.778548 -1.0882912 Nov 5 5536.755 134.311708 -1.0949055 Dec 5 5515.148 112.960656 -1.0879539 Jan 6 5492.962 94.537532 0.9175360 Feb 6 5694.161 109.042743 0.6689765 Mar 6 5849.748 115.417099 0.6624741 Apr 6 6115.996 136.073162 0.6442822 May 6 6174.348 125.429675 0.6523725 Jun 6 6512.947 154.621461 0.6332215 Jul 6 6383.125 115.669790 0.6552763 Aug 6 6673.016 139.527524 0.6436176 Sep 6 6935.974 156.429744 0.6364889 Oct 6 7300.054 184.864961 0.6261383 Nov 6 7392.300 172.181938 0.6301228 Dec 6 7496.677 162.896953 0.6326404 Jan 7 7595.989 154.222184 -11.2785124 Feb 7 7163.299 74.325154 -2.5085800 Mar 7 7198.694 68.993634 -2.5039174 Apr 7 7248.132 66.315543 -2.5019007 May 7 7350.015 71.186308 -2.5050664 Jun 7 7428.256 72.152326 -2.5056082 Jul 7 7781.034 110.581288 -2.5242126 Aug 7 7824.850 101.438538 -2.5203925 Sep 7 8183.753 136.695291 -2.5331066 Oct 7 8374.005 144.029289 -2.5353892 Nov 7 8350.239 121.051701 -2.5292171 Dec 7 8674.646 148.898713 -2.5356729 Jan 8 8777.569 142.623016 25.2211498 Feb 8 9139.600 172.514009 -1.1404687 Mar 8 9124.411 146.807886 -1.1208802 Apr 8 9024.309 112.994451 -1.0986495 May 8 8851.486 73.853751 -1.0764399 Jun 8 8865.652 65.680028 -1.0724370 Jul 8 9164.826 97.654640 -1.0859518 Aug 8 8517.709 -4.333840 -1.0487470 Sep 8 8554.491 1.296480 -1.0505196 Oct 8 7556.213 -135.583617 -1.0133256 Nov 8 7852.247 -76.478618 -1.0271867 Dec 8 7443.018 -122.044951 -1.0179639 Jan 9 7949.433 -36.224652 43.0965082 Feb 9 8266.313 11.908543 -2.2733230 Mar 9 7519.593 -91.985262 -2.2031112 Apr 9 7202.585 -122.801149 -2.1851344 May 9 7195.883 -106.902177 -2.1931393 Jun 9 6195.720 -229.226475 -2.1399841 Jul 9 5106.306 -347.020669 -2.0958066 Aug 9 4802.558 -341.094949 -2.0977247 Sep 9 4463.708 -340.787516 -2.0978106 Oct 9 4400.697 -302.749137 -2.1069818 Nov 9 4245.741 -282.510521 -2.1111932 Dec 9 4295.939 -236.949848 -2.1193756 > m$resid Jan Feb Mar Apr May Jun 1 0.00000000 0.47414318 -0.07348247 0.92151011 0.19210747 0.02864014 2 0.54484930 0.24284408 0.47727333 -0.86527639 -0.32490683 -0.21743655 3 0.28973107 -0.52916068 -0.91024314 -1.09190368 -0.30270512 -0.10543526 4 1.11018306 0.58239673 0.97645256 0.15748162 0.37925767 0.05340968 5 0.31578660 -0.80039301 0.37386594 -0.33471018 0.16838060 0.33095955 6 -0.51857793 0.37844433 0.17236841 0.55859560 -0.28784266 0.78949201 7 -0.24258719 -2.09615700 -0.14417453 -0.07242481 0.13172779 0.02612635 8 -0.17467627 0.78741058 -0.69516240 -0.91445080 -1.05856011 -0.22106435 9 2.38027621 1.27191057 -2.80962773 -0.83339855 0.42999284 -3.30838015 Jul Aug Sep Oct Nov Dec 1 -1.26741587 -0.56796473 -0.41529261 0.18143146 0.15645245 -1.19624487 2 -1.43345910 0.93703382 1.76431045 0.07032411 0.30579783 -0.05734589 3 -0.87361741 1.42620230 0.80344219 0.52063691 -1.01353613 0.42852770 4 -0.10817343 -0.08109884 0.14920641 0.79151873 0.10512391 -0.60305787 5 0.11021527 0.92539437 0.20384814 0.21737898 0.47422111 -0.57748828 6 -1.05347679 0.64526315 0.45714923 0.76908732 -0.34304059 -0.25113467 7 1.03935167 -0.24727957 0.95358476 0.19836357 -0.62148222 0.75319112 8 0.86479443 -2.75845109 0.15228325 -3.70222672 1.59863618 -1.23245629 9 -3.18591650 0.16027194 0.00831516 1.02883541 0.54740296 1.23230508 > 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/1snot1259863072.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/2ayfj1259863072.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/3o4rm1259863072.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/4fbp51259863072.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/5rjd71259863072.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/69h621259863072.tab") > system("convert tmp/1snot1259863072.ps tmp/1snot1259863072.png") > system("convert tmp/2ayfj1259863072.ps tmp/2ayfj1259863072.png") > system("convert tmp/3o4rm1259863072.ps tmp/3o4rm1259863072.png") > system("convert tmp/4fbp51259863072.ps tmp/4fbp51259863072.png") > system("convert tmp/5rjd71259863072.ps tmp/5rjd71259863072.png") > > > proc.time() user system elapsed 1.828 0.831 2.244