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(267413,267366,264777,258863,254844,254868,277267,285351,286602,283042,276687,277915,277128,277103,275037,270150,267140,264993,287259,291186,292300,288186,281477,282656,280190,280408,276836,275216,274352,271311,289802,290726,292300,278506,269826,265861,269034,264176,255198,253353,246057,235372,258556,260993,254663,250643,243422,247105,248541,245039,237080,237085,225554,226839,247934,248333,246969,245098,246263,255765,264319,268347,273046,273963,267430,271993,292710) > 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 10384851 4047765 0 > m$fitted level slope sea Jan 1 267413.0 0.000000 0.000000 Feb 1 267384.1 -30.375380 -18.122333 Mar 1 265472.6 -1235.301307 -695.599309 Apr 1 260362.3 -3831.511382 -1499.338886 May 1 254957.8 -4881.918488 -113.818656 Jun 1 253458.4 -2664.063490 1409.577840 Jul 1 270325.6 10177.982855 6941.371196 Aug 1 285820.2 13682.097697 -469.199420 Sep 1 290713.0 7887.915883 -4111.019426 Oct 1 286637.5 2.822738 -3595.495075 Nov 1 278452.6 -5392.456042 -1765.649227 Dec 1 275928.1 -3502.788864 1986.884691 Jan 2 275910.7 -1214.396079 1217.316600 Feb 2 276148.1 -256.840860 954.892717 Mar 2 274331.3 -1258.143972 705.742173 Apr 2 270423.9 -2934.228221 -273.947923 May 2 267842.6 -2708.113882 -702.636211 Jun 2 269591.2 129.630269 -4598.178829 Jul 2 279833.3 6537.525350 7425.673442 Aug 2 289111.7 8277.311256 2074.295001 Sep 2 293779.4 5982.932018 -1479.439462 Oct 2 291657.2 830.263640 -3471.212801 Nov 2 286081.2 -3240.235700 -4604.180237 Dec 2 281914.5 -3828.206167 741.541801 Jan 3 279035.0 -3225.938856 1155.005188 Feb 3 277669.2 -2043.518146 2738.828719 Mar 3 274702.6 -2625.487591 2133.362316 Apr 3 274000.7 -1419.592801 1215.285202 May 3 275819.2 620.622798 -1467.193210 Jun 3 279435.2 2508.346802 -8124.244399 Jul 3 283706.5 3615.318555 6095.482107 Aug 3 287664.3 3830.252966 3061.742171 Sep 3 290980.7 3507.362554 1319.294708 Oct 3 283684.6 -3284.537960 -5178.568397 Nov 3 275220.7 -6539.310947 -5394.719885 Dec 3 266141.3 -8134.963078 -280.309029 Jan 4 265229.0 -3594.223929 3805.048474 Feb 4 261316.3 -3794.396000 2859.737898 Mar 4 254440.7 -5723.069614 757.267352 Apr 4 251795.2 -3801.651170 1557.797376 May 4 248835.1 -3274.980990 -2778.148362 Jun 4 245060.1 -3588.276266 -9688.110371 Jul 4 250197.8 1868.690572 8358.157479 Aug 4 256704.0 4765.631498 4289.011771 Sep 4 252651.5 -744.879260 2011.469902 Oct 4 252537.9 -350.253529 -1894.916556 Nov 4 249500.1 -2030.277331 -6078.056650 Dec 4 249297.7 -887.439653 -2192.716621 Jan 5 245541.8 -2681.882525 2999.234435 Feb 5 241233.4 -3698.700729 3805.575664 Mar 5 237078.6 -3983.132145 1.393087 Apr 5 234613.1 -3038.061132 2471.865894 May 5 229446.7 -4365.186192 -3892.719834 Jun 5 236093.5 2508.794200 -9254.543571 Jul 5 241002.6 4005.780358 6931.367944 Aug 5 242033.6 2152.655578 6299.406979 Sep 5 244729.4 2490.923308 2239.587358 Oct 5 246695.1 2163.780594 -1597.128901 Nov 5 251948.3 4088.505533 -5685.266906 Dec 5 256798.3 4563.142932 -1033.299902 Jan 6 260745.2 4178.958525 3573.754715 Feb 6 264211.3 3734.810420 4135.695924 Mar 6 271491.4 5939.701333 1554.619764 Apr 6 273021.2 3199.465454 941.791858 May 6 275571.6 2795.775298 -8141.588732 Jun 6 281004.5 4437.433766 -9011.530487 Jul 6 284960.1 4137.619062 7749.874722 > m$resid Jan Feb Mar Apr May Jun 1 0.00000000 -0.01186856 -0.38180737 -0.85961113 -0.32291222 0.68682435 2 0.71284573 0.30202125 -0.30801614 -0.52764020 0.07063611 0.87683729 3 0.18749256 0.36755088 -0.18004148 0.37558580 0.63612313 0.58508378 4 1.41136746 -0.06211187 -0.59761165 0.59716201 0.16389373 -0.09722729 5 -0.55726086 -0.31539133 -0.08818646 0.29352703 -0.41258498 2.13409355 6 -0.11925577 -0.13775249 0.68380694 -0.85085125 -0.12543888 0.50971815 Jul Aug Sep Oct Nov Dec 1 3.99268796 1.08745140 -1.79792560 -2.44692621 -1.67422917 0.58638743 2 1.98777506 0.54033575 -0.71196308 -1.59896466 -1.26318802 -0.18254672 3 0.34300676 0.06670888 -0.10021273 -2.10761970 -1.01010820 -0.49555403 4 1.69138685 0.89861411 -1.71012709 0.12247373 -0.52150058 0.35491313 5 0.46418648 -0.57471500 0.10496350 -0.10154299 0.59757264 0.14738679 6 -0.09299350 > 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/18czu1259963910.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/2jmlr1259963910.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/32vhl1259963910.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/4ydp51259963910.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/51qql1259963910.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/6a2p91259963910.tab") > system("convert tmp/18czu1259963910.ps tmp/18czu1259963910.png") > system("convert tmp/2jmlr1259963910.ps tmp/2jmlr1259963910.png") > system("convert tmp/32vhl1259963910.ps tmp/32vhl1259963910.png") > system("convert tmp/4ydp51259963910.ps tmp/4ydp51259963910.png") > system("convert tmp/51qql1259963910.ps tmp/51qql1259963910.png") > > > proc.time() user system elapsed 1.300 0.820 2.969