R version 2.12.0 (2010-10-15) Copyright (C) 2010 The R Foundation for Statistical Computing ISBN 3-900051-07-0 Platform: i486-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(989236,1008380,1207763,1368839,1469798,1498721,1761769,1653214,1599104,1421179,1163995,1037735,1015407,1039210,1258049,1469445,1552346,1549144,1785895,1662335,1629440,1467430,1202209,1076982,1039367,1063449,1335135,1491602,1591972,1641248,1898849,1798580,1762444,1622044,1368955,1262973,1195650,1269530,1479279,1607819,1712466,1721766,1949843,1821326,1757802,1590367,1260647,1149235,1016367,1027885,1262159,1520854,1544144,1564709,1821776,1741365,1623386,1498658,1241822,1136029) > 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 8342732573 278650223 0 > m$fitted level slope sea Jan 1 989236 0.000 0.0000 Feb 1 1007465 18304.458 915.1760 Mar 1 1190887 164849.202 16875.8458 Apr 1 1371097 178472.562 -2257.5587 May 1 1475942 113401.039 -6144.2596 Jun 1 1503787 37693.489 -5065.8174 Jul 1 1741128 214372.611 20641.2750 Aug 1 1681601 -27997.941 -28387.1570 Sep 1 1596785 -78275.316 2319.4886 Oct 1 1427149 -159118.312 -5969.8137 Nov 1 1169752 -246083.611 -5756.5881 Dec 1 1025503 -155971.650 12232.0308 Jan 2 1006999 -34428.922 8408.2219 Feb 2 1050488 34586.329 -11278.3930 Mar 2 1233258 160152.959 24791.1547 Apr 2 1459806 216956.596 9638.9626 May 2 1555364 113678.637 -3018.2469 Jun 2 1596435 51863.311 -47291.2106 Jul 2 1725617 117714.651 60278.1890 Aug 2 1706520 1195.996 -44185.0110 Sep 2 1628079 -66624.214 1360.5167 Oct 2 1465089 -148692.705 2341.4063 Nov 2 1221073 -229874.489 -18864.4675 Dec 2 1064616 -167362.510 12366.4465 Jan 3 1020007 -62890.791 19359.8632 Feb 3 1083918 45164.879 -20468.5366 Mar 3 1309673 196678.198 25461.6181 Apr 3 1472032 167735.585 19570.2414 May 3 1588490 124610.668 3481.8545 Jun 3 1707722 120088.898 -66473.9333 Jul 3 1823943 116835.559 74906.0112 Aug 3 1848905 39552.856 -50324.6589 Sep 3 1764562 -64662.477 -2118.0813 Oct 3 1609138 -141007.301 12906.3489 Nov 3 1396425 -201320.553 -27469.8645 Dec 3 1250238 -154962.615 12734.5804 Jan 4 1178212 -85223.319 17438.3381 Feb 4 1299014 88010.057 -29483.7494 Mar 4 1446328 137524.819 32951.4316 Apr 4 1581536 135585.953 26282.6535 May 4 1711512 130893.193 954.2228 Jun 4 1799807 95297.588 -78040.9936 Jul 4 1871248 75353.376 78594.6805 Aug 4 1867720 9396.032 -46393.9215 Sep 4 1761259 -87478.384 -3457.2365 Oct 4 1572182 -172429.180 18184.5908 Nov 4 1297721 -257731.673 -37074.3850 Dec 4 1126799 -185181.565 22436.4014 Jan 5 1019681 -119917.972 -3313.9017 Feb 5 1050169 5746.213 -22284.3079 Mar 5 1213281 136714.031 48878.4319 Apr 5 1481909 246666.669 38945.1247 May 5 1559352 105621.747 -15207.8525 Jun 5 1647632 91182.532 -82923.0650 Jul 5 1739556 91799.809 82220.2625 Aug 5 1781924 50619.564 -40558.5752 Sep 5 1634953 -113991.417 -11567.4194 Oct 5 1463630 -161752.955 35028.1310 Nov 5 1286067 -174921.488 -44244.6498 Dec 5 1113060 -173326.691 22968.6152 > m$resid Jan Feb Mar Apr May Jun 1 0.000000000 0.205158157 1.683816167 0.148470105 -0.712754757 -0.828882468 2 1.332186563 0.763051784 1.389874132 0.622455116 -1.129893240 -0.676935426 3 1.145478655 1.185080709 1.663346354 -0.317588360 -0.471632508 -0.049514801 4 0.764564637 1.896581827 0.542615429 -0.021268445 -0.051335702 -0.389706535 5 0.715242110 1.375331445 1.434508436 1.205546234 -1.543425779 -0.158059884 Jul Aug Sep Oct Nov Dec 1 1.934336890 -2.653541526 -0.550450694 -0.885091652 -0.952120316 0.986570853 2 0.720946036 -1.275683946 -0.742514535 -0.898509046 -0.888806448 0.684409798 3 -0.035619190 -0.846106048 -1.140980314 -0.835848920 -0.660328170 0.507591916 4 -0.218370851 -0.722108633 -1.060611369 -0.930079379 -0.933906473 0.794472369 5 0.006758736 -0.450848753 -1.802212860 -0.522916606 -0.144171436 0.017465305 > 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 > postscript(file="/var/www/rcomp/tmp/1sugm1292604146.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') > grid() > dev.off() null device 1 > mylagmax <- nx/2 > postscript(file="/var/www/rcomp/tmp/2k3f71292604146.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/rcomp/tmp/3k3f71292604146.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/rcomp/tmp/40y3g1292604146.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/rcomp/tmp/564vv1292604146.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/rcomp/createtable file can be downloaded at http://www.wessa.net/cretab > load(file="/var/www/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/rcomp/tmp/694u01292604146.tab") > > try(system("convert tmp/1sugm1292604146.ps tmp/1sugm1292604146.png",intern=TRUE)) character(0) > try(system("convert tmp/2k3f71292604146.ps tmp/2k3f71292604146.png",intern=TRUE)) character(0) > try(system("convert tmp/3k3f71292604146.ps tmp/3k3f71292604146.png",intern=TRUE)) character(0) > try(system("convert tmp/40y3g1292604146.ps tmp/40y3g1292604146.png",intern=TRUE)) character(0) > try(system("convert tmp/564vv1292604146.ps tmp/564vv1292604146.png",intern=TRUE)) character(0) > > > proc.time() user system elapsed 1.500 0.490 1.978