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(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) > 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 9557919.081 1583.144 2694758.760 0.000 > m$fitted level slope sea Jan 1 17409.00 0.0000 0.0000 Feb 1 12506.29 -244.3877 -992.2913 Mar 1 25092.05 263.4559 6421.9516 Apr 1 27652.21 302.4626 -581.2129 May 1 29050.32 311.5627 411.6819 Jun 1 27110.52 295.1162 -1005.5237 Jul 1 23562.94 264.3724 -1165.9421 Aug 1 23171.24 258.7947 671.7647 Sep 1 22014.81 246.4822 -309.8068 Oct 1 19115.24 218.8415 -1026.2362 Nov 1 19741.62 222.4514 1022.3830 Dec 1 23451.97 253.6226 1864.0325 Jan 2 20987.76 339.9959 -3283.7645 Feb 2 21350.34 340.2216 -5802.3356 Mar 2 22150.69 350.7073 5878.3084 Apr 2 27461.34 439.1513 1921.6570 May 2 32981.00 495.5257 3457.0038 Jun 2 32995.20 491.7072 -961.2005 Jul 2 27228.82 445.7067 -4549.8153 Aug 2 23951.94 416.9424 367.0640 Sep 2 19245.39 374.5238 -1241.3871 Oct 2 18338.76 363.5486 -801.7647 Nov 2 19417.43 368.8802 948.5658 Dec 2 19857.96 369.0774 2924.0447 Jan 3 21581.04 362.7887 -2412.0356 Feb 3 20963.93 356.2291 -7156.9300 Mar 3 24394.84 404.3695 5348.1560 Apr 3 25753.52 419.7107 -162.5234 May 3 25789.80 414.8295 3306.1963 Jun 3 25505.02 407.8874 976.9797 Jul 3 25465.32 403.9746 -3060.3241 Aug 3 25071.31 397.0657 1972.6881 Sep 3 21178.02 358.9875 -3208.0215 Oct 3 19966.76 345.4938 -1236.7605 Nov 3 19136.39 337.2463 547.6083 Dec 3 17904.53 330.7069 1880.4685 Jan 4 19608.73 334.6352 -1129.7327 Feb 4 19524.10 331.3918 -8826.0990 Mar 4 24308.10 389.2641 7647.9031 Apr 4 28166.33 439.6313 1339.6682 May 4 30461.58 464.1203 4044.4203 Jun 4 28128.72 432.4665 -963.7239 Jul 4 28755.73 434.4353 -2019.7281 Aug 4 23498.52 379.4943 192.4796 Sep 4 21292.41 355.1571 -3135.4082 Oct 4 18991.86 331.7573 -1663.8576 Nov 4 17557.77 318.4340 647.2313 Dec 4 18658.80 323.2443 2336.1973 Jan 5 18869.78 322.5433 -1487.7802 Feb 5 19849.13 328.4289 -10482.1267 > m$resid Jan Feb Mar Apr May Jun 1 0.00000000 -1.02575577 3.76129702 0.74187558 0.35536483 -0.72944554 2 -0.97667465 0.00699642 0.13860059 1.57609702 1.64104500 -0.15562709 3 0.44832968 -0.31240153 0.95332870 0.30184415 -0.12323146 -0.22570751 4 0.44589012 -0.13394527 1.39988120 1.09826822 0.59430512 -0.90029561 5 -0.03618418 0.20980169 Jul Aug Sep Oct Nov Dec 1 -1.24419096 -0.21235623 -0.45800566 -1.01803451 0.13185972 1.12836445 2 -2.02217824 -1.20227708 -1.65448116 -0.41376364 0.23080378 0.02313670 3 -0.14441303 -0.25740047 -1.38367004 -0.50618790 -0.37865637 -0.50631257 4 0.06267970 -1.83369999 -0.83268998 -0.85453856 -0.56779657 0.25201693 5 > 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/15p2c1259785002.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/2x7lo1259785002.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/3tjhi1259785002.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/43ac71259785002.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/51giy1259785002.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/63mp71259785002.tab") > system("convert tmp/15p2c1259785002.ps tmp/15p2c1259785002.png") > system("convert tmp/2x7lo1259785002.ps tmp/2x7lo1259785002.png") > system("convert tmp/3tjhi1259785002.ps tmp/3tjhi1259785002.png") > system("convert tmp/43ac71259785002.ps tmp/43ac71259785002.png") > system("convert tmp/51giy1259785002.ps tmp/51giy1259785002.png") > > > proc.time() user system elapsed 1.336 0.808 1.590