R version 2.15.2 (2012-10-26) -- "Trick or Treat" Copyright (C) 2012 The R Foundation for Statistical Computing ISBN 3-900051-07-0 Platform: i686-pc-linux-gnu (32-bit) 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' > 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 54604452029 15942638858 0 > m$fitted level slope sea Jan 1 655362.0 0.000 0.000 Feb 1 804406.3 154723.670 68720.729 Mar 1 1106300.6 252402.524 1596.362 Apr 1 1514241.7 361658.416 41722.265 May 1 1734860.3 263973.661 -63701.348 Jun 1 1606381.3 -4624.909 -113073.341 Jul 1 2589589.2 673905.401 368206.769 Aug 1 2870691.9 403718.870 -232000.933 Sep 1 1764653.9 -634755.496 -458984.853 Oct 1 1125357.5 -637878.216 155138.514 Nov 1 813555.9 -413672.625 108344.120 Dec 1 776975.1 -154394.443 90912.852 Jan 2 674263.0 -118967.082 -21677.043 Feb 2 817118.3 61258.582 96712.703 Mar 2 1100837.2 209625.373 7706.803 Apr 2 1402058.8 270175.012 153768.215 May 2 1685824.3 279252.433 13458.692 Jun 2 2025952.9 319595.259 -516494.938 Jul 2 2693984.9 549922.935 574990.074 Aug 2 2455590.5 27803.918 -30574.536 Sep 2 1887240.7 -367364.789 -574537.742 Oct 2 1322291.5 -498337.762 43206.512 Nov 2 883149.6 -459116.921 51303.426 Dec 2 614397.0 -333135.170 160621.953 Jan 3 633505.4 -99975.501 17636.643 Feb 3 752375.1 45074.412 90816.882 Mar 3 1071980.1 225370.487 74785.858 Apr 3 1471062.7 338941.273 181538.334 May 3 1626318.2 218206.615 -160412.199 Jun 3 2205887.1 455393.024 -553153.051 Jul 3 2258299.2 191716.792 664034.823 Aug 3 2485325.2 214831.174 217479.812 Sep 3 2139065.6 -152861.885 -680109.602 Oct 3 1487508.6 -479723.845 -77145.586 Nov 3 985788.8 -494133.655 33490.174 Dec 3 774886.6 -308681.244 161687.429 Jan 4 691302.9 -161180.411 17614.095 Feb 4 800990.6 16266.165 84304.352 Mar 4 1017114.0 146601.339 82548.995 Apr 4 1322725.0 250074.738 253495.024 May 4 1748618.9 364798.302 -260748.930 Jun 4 2010215.7 297435.388 -521580.705 Jul 4 2267215.0 271094.349 615314.961 Aug 4 2324713.9 132015.681 352312.132 Sep 4 2019831.9 -152644.417 -615433.928 Oct 4 1491863.1 -397249.207 -147493.084 Nov 4 991231.0 -464609.892 -54366.037 Dec 4 704400.3 -348773.930 168304.660 Jan 5 609708.9 -183123.028 18442.106 Feb 5 811251.9 67500.559 142460.098 Mar 5 1091172.4 205525.394 69211.629 Apr 5 1248644.2 174338.291 151973.825 May 5 1797784.6 417973.358 -136273.559 Jun 5 2109124.0 348610.229 -613777.001 Jul 5 2296705.4 243985.537 622080.625 Aug 5 2327785.8 105763.616 447891.211 Sep 5 1999944.6 -175778.823 -592918.584 Oct 5 1527005.3 -368766.836 -156806.294 Nov 5 1068224.4 -427227.539 -103698.383 Dec 5 726189.3 -371883.920 124661.726 Jan 6 681906.5 -158988.664 1211.505 Feb 6 711405.3 -36586.682 135818.698 Mar 6 929423.3 128451.266 143832.657 Apr 6 1412776.7 358337.190 101549.295 May 6 1684675.4 302288.604 -180941.447 Jun 6 2069108.5 355589.320 -561396.477 Jul 6 2249922.9 242260.613 615775.097 Aug 6 2252724.5 87120.402 535403.471 Sep 6 1969764.6 -152592.210 -578168.571 Oct 6 1535961.7 -334761.550 -169583.710 Nov 6 1084296.5 -410508.113 -138001.505 Dec 6 782388.2 -340116.398 77237.760 > m$resid Jan Feb Mar Apr May Jun 1 0.00000000 0.79269778 0.43495734 0.49012522 -0.41379765 -1.14966331 2 0.15213384 0.78369778 0.63019621 0.26321687 0.03896431 0.17202078 3 1.00094352 0.62171716 0.76940526 0.48828736 -0.51853430 1.01320209 4 0.63231491 0.75930926 0.55701014 0.44369734 0.49214338 -0.28812488 5 0.70949258 1.07209031 0.59021366 -0.13361239 1.04435181 -0.29684398 6 0.91142083 0.52356638 0.70591612 0.98450814 -0.24014974 0.22815425 Jul Aug Sep Oct Nov Dec 1 2.90713343 -1.15615649 -4.44410989 -0.01336369 0.95947078 1.10956307 2 0.98611830 -2.23553728 -1.69095076 -0.56052322 0.16784762 0.53935499 3 -1.12747130 0.09895119 -1.57359275 -1.39878280 -0.06166914 0.79421386 4 -0.11261605 -0.59513341 -1.21828653 -1.04681171 -0.28832010 0.49606459 5 -0.44739591 -0.59132096 -1.20487952 -0.82599799 -0.25026225 0.23698280 6 -0.48472218 -0.66363497 -1.02579738 -0.77975995 -0.32429438 0.30139084 > 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/fisher/rcomp/tmp/1q5u21354484395.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/fisher/rcomp/tmp/2g5fu1354484396.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/fisher/rcomp/tmp/3mjfg1354484396.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/fisher/rcomp/tmp/4xuko1354484396.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/fisher/rcomp/tmp/5o1b41354484396.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/fisher/rcomp/createtable file can be downloaded at http://www.wessa.net/cretab > load(file="/var/fisher/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/fisher/rcomp/tmp/6vzwt1354484396.tab") > > try(system("convert tmp/1q5u21354484395.ps tmp/1q5u21354484395.png",intern=TRUE)) character(0) > try(system("convert tmp/2g5fu1354484396.ps tmp/2g5fu1354484396.png",intern=TRUE)) character(0) > try(system("convert tmp/3mjfg1354484396.ps tmp/3mjfg1354484396.png",intern=TRUE)) character(0) > try(system("convert tmp/4xuko1354484396.ps tmp/4xuko1354484396.png",intern=TRUE)) character(0) > try(system("convert tmp/5o1b41354484396.ps tmp/5o1b41354484396.png",intern=TRUE)) character(0) > > > proc.time() user system elapsed 2.950 0.713 3.666