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(655362,873127,1107897,1555964,1671159,1493308,2957796,2638691,1305669,1280496,921900,867888,652586,913831,1108544,1555827,1699283,1509458,3268975,2425016,1312703,1365498,934453,775019,651142,843192,1146766,1652601,1465906,1652734,2922334,2702805,1458956,1410363,1019279,936574,708917,885295,1099663,1576220,1487870,1488635,2882530,2677026,1404398,1344370,936865,872705,628151,953712,1160384,1400618,1661511,1495347,2918786,2775677,1407026,1370199,964526,850851,683118,847224,1073256,1514326,1503734,1507712,2865698,2788128,1391596,1366378,946295,859626) > 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 54604455798 15942635265 0 > m$fitted level slope sea Jan 1 655362.0 0.000 0.000 Feb 1 804406.3 154723.683 68720.714 Mar 1 1106300.6 252402.524 1596.359 Apr 1 1514241.7 361658.422 41722.261 May 1 1734860.3 263973.653 -63701.343 Jun 1 1606381.3 -4624.931 -113073.325 Jul 1 2589589.3 673905.456 368206.737 Aug 1 2870691.9 403718.853 -232000.932 Sep 1 1764653.8 -634755.587 -458984.791 Oct 1 1125357.5 -637878.217 155138.535 Nov 1 813555.9 -413672.590 108344.104 Dec 1 776975.2 -154394.415 90912.826 Jan 2 674263.1 -118967.083 -21677.052 Feb 2 817118.3 61258.591 96712.693 Mar 2 1100837.2 209625.384 7706.787 Apr 2 1402058.8 270175.015 153768.205 May 2 1685824.3 279252.418 13458.707 Jun 2 2025952.9 319595.248 -516494.912 Jul 2 2693985.0 549922.976 574990.047 Aug 2 2455590.5 27803.900 -30574.547 Sep 2 1887240.7 -367364.841 -574537.710 Oct 2 1322291.5 -498337.788 43206.550 Nov 2 883149.6 -459116.910 51303.437 Dec 2 614397.1 -333135.144 160621.938 Jan 3 633505.4 -99975.477 17636.619 Feb 3 752375.1 45074.424 90816.862 Mar 3 1071980.2 225370.496 74785.841 Apr 3 1471062.7 338941.270 181538.333 May 3 1626318.2 218206.590 -160412.171 Jun 3 2205887.0 455393.041 -553153.048 Jul 3 2258299.2 191716.798 664034.809 Aug 3 2485325.2 214831.186 217479.785 Sep 3 2139065.6 -152861.924 -680109.586 Oct 3 1487508.5 -479723.888 -77145.546 Nov 3 985788.8 -494133.659 33490.198 Dec 3 774886.6 -308681.220 161687.424 Jan 4 691302.9 -161180.388 17614.077 Feb 4 800990.7 16266.182 84304.330 Mar 4 1017114.0 146601.345 82548.981 Apr 4 1322725.0 250074.735 253495.025 May 4 1748618.9 364798.299 -260748.917 Jun 4 2010215.7 297435.388 -521580.699 Jul 4 2267215.1 271094.361 615314.947 Aug 4 2324713.9 132015.685 352312.109 Sep 4 2019831.9 -152644.443 -615433.926 Oct 4 1491863.1 -397249.244 -147493.053 Nov 4 991231.0 -464609.908 -54366.007 Dec 4 704400.3 -348773.916 168304.667 Jan 5 609708.9 -183123.004 18442.095 Feb 5 811251.9 67500.584 142460.076 Mar 5 1091172.4 205525.402 69211.618 Apr 5 1248644.2 174338.279 151973.834 May 5 1797784.6 417973.371 -136273.560 Jun 5 2109124.0 348610.229 -613777.000 Jul 5 2296705.4 243985.542 622080.612 Aug 5 2327785.8 105763.618 447891.189 Sep 5 1999944.6 -175778.847 -592918.586 Oct 5 1527005.3 -368766.867 -156806.273 Nov 5 1068224.4 -427227.558 -103698.355 Dec 5 726189.3 -371883.918 124661.744 Jan 6 681906.5 -158988.640 1211.502 Feb 6 711405.3 -36586.664 135818.688 Mar 6 929423.4 128451.278 143832.649 Apr 6 1412776.7 358337.200 101549.292 May 6 1684675.4 302288.601 -180941.442 Jun 6 2069108.5 355589.330 -561396.486 Jul 6 2249922.9 242260.618 615775.078 Aug 6 2252724.6 87120.399 535403.448 Sep 6 1969764.6 -152592.232 -578168.577 Oct 6 1535961.7 -334761.579 -169583.694 Nov 6 1084296.5 -410508.134 -138001.477 Dec 6 782388.2 -340116.398 77237.782 > m$resid Jan Feb Mar Apr May Jun 1 0.00000000 0.79269779 0.43495727 0.49012522 -0.41379770 -1.14966333 2 0.15213371 0.78369780 0.63019621 0.26321682 0.03896423 0.17202079 3 1.00094348 0.62171709 0.76940523 0.48828729 -0.51853437 1.01320224 4 0.63231488 0.75930921 0.55701008 0.44369728 0.49214337 -0.28812485 5 0.70949260 1.07209028 0.59021356 -0.13361247 1.04435188 -0.29684402 6 0.91142088 0.52356633 0.70591608 0.98450810 -0.24014979 0.22815430 Jul Aug Sep Oct Nov Dec 1 2.90713366 -1.15615676 -4.44411005 -0.01336330 0.95947091 1.10956300 2 0.98611849 -2.23553746 -1.69095084 -0.56052309 0.16784778 0.53935504 3 -1.12747131 0.09895121 -1.57359291 -1.39878277 -0.06166897 0.79421396 4 -0.11261599 -0.59513343 -1.21828662 -1.04681172 -0.28831999 0.49606470 5 -0.44739588 -0.59132095 -1.20487959 -0.82599799 -0.25026219 0.23698288 6 -0.48472219 -0.66363497 -1.02579743 -0.77975996 -0.32429433 0.30139091 > 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/1fsdf1293532966.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/html/rcomp/tmp/28kci1293532966.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/html/rcomp/tmp/3jtt31293532966.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/html/rcomp/tmp/4jtt31293532966.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/html/rcomp/tmp/5tks61293532966.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/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/6flru1293532966.tab") > > try(system("convert tmp/1fsdf1293532966.ps tmp/1fsdf1293532966.png",intern=TRUE)) character(0) > try(system("convert tmp/28kci1293532966.ps tmp/28kci1293532966.png",intern=TRUE)) character(0) > try(system("convert tmp/3jtt31293532966.ps tmp/3jtt31293532966.png",intern=TRUE)) character(0) > try(system("convert tmp/4jtt31293532966.ps tmp/4jtt31293532966.png",intern=TRUE)) character(0) > try(system("convert tmp/5tks61293532966.ps tmp/5tks61293532966.png",intern=TRUE)) character(0) > > > proc.time() user system elapsed 1.347 0.816 3.033