R version 2.13.0 (2011-04-13) Copyright (C) 2011 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(127,128.7,127.3,136.7,133.8,137.2,147.4,137.6,123.6,117.4,113.7,106.8,103.3,96.3,96.2,94.7,94.6,95.7,106.7,100.2,94.2,97.6,94.3,98,93.6,86.3,90.7,81.8,87.6,73.8,63.3,59.1,52.9,54.6,52.4,67.5,90.4,126,144.3,167.8,166.2,156,137,129.3,118,114.7,112.8,115.7,103.9,96.9,88.8,93,86.3,82.3,82.4,76.6,72.7,67.5,77.3,73.7,73,78.2,90.7,91.5,86.3,86.8,86.1,77.1,75.7,78.7,71.5,69.6,73.6,78.1,78.3,71.5,68.7,61.2,64.7,64.6,56.3,54.5,49.5,54,59.2,52.4,52.8,47.8,45.2,47.1,42.6,42.1,39.4,39.6,37.8,36.8,36.5,34,37.5,36.1,35.1,32.8,32.2,38.1,41.4,38.5,34.6,31.2,34.7,35.8,33.8,32.5,31.2,32.5,32.3,30.1,25.3,24.5,23.7,23.8,26.9,32.4,32.6,31.7,34.1,34.9,33.1,32.3,34.6,32.7,32.6,41,40) > 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 21.05990 14.56139 0.00000 0.00000 > m$fitted level slope sea Jan 1 127.00000 0.00000000 0.0000000000 Feb 1 128.62778 0.38335753 0.0722208948 Mar 1 127.20848 -0.59051936 0.0915240857 Apr 1 136.65681 4.93500178 0.0431946236 May 1 133.73992 0.59074989 0.0600812396 Jun 1 137.14261 2.14916214 0.0573853348 Jul 1 147.34605 6.61574959 0.0539456989 Aug 1 137.54293 -2.49185076 0.0570680289 Sep 1 123.54196 -8.87669657 0.0580424982 Oct 1 117.34206 -7.39163519 0.0579415948 Nov 1 113.64212 -5.34347413 0.0578796407 Dec 1 106.74211 -6.20705236 0.0578912700 Jan 2 103.15670 -4.79183370 0.1432965530 Feb 2 96.35447 -5.77713415 -0.0544677470 Mar 2 96.27684 -2.59060349 -0.0768397819 Apr 2 94.77875 -1.98412750 -0.0787475968 May 2 94.68021 -0.93831418 -0.0802140769 Jun 2 95.78092 0.19266816 -0.0809202839 Jul 2 106.78259 6.18880630 -0.0825872000 Aug 2 100.28172 -0.85111262 -0.0817159192 Sep 2 94.28156 -3.70777584 -0.0815585229 Oct 2 97.68166 0.23569948 -0.0816552523 Nov 2 94.38163 -1.72594681 -0.0816338311 Dec 2 98.08165 1.28443013 -0.0816484660 Jan 3 93.51285 -1.92214235 0.0871525868 Feb 3 86.37671 -4.59899728 -0.0767079655 Mar 3 90.79969 0.43154084 -0.0996884849 Apr 3 81.88910 -4.75345646 -0.0890995572 May 3 87.69443 1.10345809 -0.0944289389 Jun 3 73.89108 -7.16563588 -0.0910786566 Jul 3 63.39074 -9.01560052 -0.0907449651 Aug 3 59.19096 -6.34384448 -0.0909595135 Sep 3 52.99096 -6.26403789 -0.0909623666 Oct 3 54.69103 -1.84549881 -0.0910326894 Nov 3 52.49103 -2.04218010 -0.0910312958 Dec 3 67.59106 7.46849617 -0.0910612955 Jan 4 88.21077 14.70509732 2.1892267281 Feb 4 125.97534 26.79656963 0.0246613623 Mar 4 144.25929 22.05618349 0.0407101749 Apr 4 167.76050 22.85811694 0.0394977488 May 4 166.15136 9.28549649 0.0486387021 Jun 4 155.94812 -1.52574871 0.0518806085 Jul 4 136.94683 -11.22074542 0.0531748847 Aug 4 129.24694 -9.26739253 0.0530587917 Sep 4 117.94691 -10.39510663 0.0530886296 Oct 4 114.64696 -6.45866061 0.0530422616 Nov 4 112.74697 -3.92946389 0.0530289987 Dec 4 115.64698 -0.14039974 0.0530201529 Jan 5 105.17517 -5.83751668 -1.2751672003 Feb 5 96.80351 -7.18246984 0.0964876629 Mar 5 88.70214 -7.69378738 0.0978627505 Apr 5 92.91006 -1.08917558 0.0899359040 May 5 86.20840 -4.20266719 0.0916003002 Jun 5 82.20843 -4.09022125 0.0915735369 Jul 5 82.30867 -1.76542299 0.0913271982 Aug 5 76.50857 -4.00385518 0.0914327917 Sep 5 72.60857 -3.94623513 0.0914315816 Oct 5 67.40856 -4.64183816 0.0914380850 Nov 5 77.20860 3.37065778 0.0914047352 Dec 5 73.60859 -0.49674234 0.0914119014 Jan 6 73.96037 -0.02824884 -0.9603672662 Feb 6 78.08635 2.19385262 0.1136539766 Mar 6 90.59915 7.93287266 0.1008525166 Apr 6 91.39521 3.97258877 0.1047932467 May 6 86.19295 -1.11720286 0.1070488879 Jun 6 86.69313 -0.21993466 0.1068718480 Jul 6 85.99310 -0.48628187 0.1068952443 Aug 6 76.99292 -5.20979471 0.1070799629 Sep 6 75.59296 -3.09607686 0.1070431638 Oct 6 78.59298 0.28609598 0.1070169500 Nov 6 71.39297 -3.86727514 0.1070312811 Dec 6 69.49297 -2.77580854 0.1070296045 Jan 7 74.43644 1.48975406 -0.8364353847 Feb 7 78.01268 2.61194940 0.0873205035 Mar 7 78.21012 1.26956577 0.0898784917 Apr 7 71.40631 -3.21030823 0.0936855497 May 7 68.60640 -2.98263648 0.0935993855 Jun 7 61.10598 -5.48898235 0.0940216966 Jul 7 64.60635 -0.50176867 0.0936475865 Aug 7 64.50636 -0.27886263 0.0936401425 Sep 7 56.20629 -4.72908160 0.0937063050 Oct 7 54.40630 -3.10399380 0.0936955490 Nov 7 49.40630 -4.15591958 0.0936986486 Dec 7 53.90631 0.64648293 0.0936923489 Jan 8 60.02416 3.67165624 -0.8241584574 Feb 8 52.37840 -2.43998635 0.0216014940 Mar 8 52.78103 -0.86000810 0.0189736508 Apr 8 47.77932 -3.15817035 0.0206778929 May 8 45.17942 -2.84845777 0.0205756144 Jun 8 47.07981 -0.21390325 0.0201882651 Jul 8 42.57966 -2.59188847 0.0203439168 Aug 8 42.07969 -1.43128401 0.0203100968 Sep 8 39.37968 -2.13518217 0.0203192283 Oct 8 39.57969 -0.83959653 0.0203117459 Nov 8 37.77969 -1.37243929 0.0203131159 Dec 8 36.77969 -1.16580574 0.0203128794 Jan 9 36.68806 -0.57162341 -0.1880579731 Feb 9 33.99177 -1.72256494 0.0082262743 Mar 9 37.49606 1.18202163 0.0039404150 Apr 9 36.09512 -0.25120351 0.0048831639 May 9 35.09500 -0.66668397 0.0050048644 Jun 9 32.79488 -1.57288341 0.0051230413 Jul 9 32.19491 -1.03311462 0.0050917037 Aug 9 38.09501 2.81345898 0.0049922831 Sep 9 41.39501 3.08339751 0.0049891770 Oct 9 38.49499 -0.23625949 0.0050061823 Nov 9 34.59499 -2.26894435 0.0050108179 Dec 9 31.19499 -2.89646687 0.0050114550 Jan 10 34.57178 0.57476672 0.1282237484 Feb 10 35.80917 0.93460569 -0.0091733016 Mar 10 33.80701 -0.69710489 -0.0070097693 Apr 10 32.50681 -1.03174274 -0.0068119962 May 10 31.20677 -1.18058768 -0.0067728242 Jun 10 32.50693 0.19569394 -0.0069340791 Jul 10 32.30692 -0.02384229 -0.0069226276 Aug 10 30.10689 -1.23120011 -0.0068945905 Sep 10 25.30687 -3.21121109 -0.0068741208 Oct 10 24.50688 -1.87344375 -0.0068802778 Nov 10 23.70688 -1.27788496 -0.0068814981 Dec 10 23.80688 -0.51341873 -0.0068821954 Jan 11 26.73393 1.39083581 0.1660727911 Feb 11 32.40056 3.71745895 -0.0005596187 Mar 11 32.59821 1.76206146 0.0017944227 Apr 11 31.69741 0.28453913 0.0025871721 May 11 34.09769 1.45831225 0.0023067419 Jun 11 34.89765 1.09306750 0.0023455916 Jul 11 33.09758 -0.51204412 0.0024215991 Aug 11 32.29758 -0.67180541 0.0024249670 Sep 11 34.59759 0.97698619 0.0024094929 Oct 11 32.69758 -0.61919847 0.0024161620 Nov 11 32.59758 -0.33114126 0.0024156262 Dec 11 40.99759 4.51299509 0.0024116149 Jan 12 40.15428 1.54769312 -0.1542773387 > m$resid Jan Feb Mar Apr May Jun 1 0.00000000 0.19058932 -0.26181484 1.45538032 -1.13959553 0.40847693 2 0.42827072 -0.23550304 0.81938629 0.15833799 0.27386111 0.29634012 3 -0.90231976 -0.66830550 1.30220644 -1.35547532 1.53411629 -2.16678089 4 1.98834054 3.06685692 -1.23102545 0.20977599 -3.55555540 -2.83298032 5 -1.54686712 -0.34390104 -0.13303279 1.72832693 -0.81568668 0.02946578 6 0.12630283 0.57097364 1.49500125 -1.03659938 -1.33351125 0.23512598 7 1.14453611 0.28929665 -0.34999324 -1.17280412 0.05965141 -0.65678357 8 0.80897631 -1.57925047 0.41220819 -0.60172256 0.08114869 0.69038391 9 0.15849121 -0.29792748 0.75817321 -0.37529546 -0.10886338 -0.23747011 10 0.92408946 0.09327390 -0.42609008 -0.08763315 -0.03900067 0.36065656 11 0.50614060 0.60374522 -0.51078178 -0.38695064 0.30755853 -0.09571315 12 -0.78714659 Jul Aug Sep Oct Nov Dec 1 1.17055421 -2.38674424 -1.67320673 0.38917331 0.53673837 -0.22630815 2 1.57129412 -1.84486105 -0.74861247 1.03342162 -0.51406637 0.78889533 3 -0.48479008 0.70015418 0.02091399 1.15791625 -0.05154203 2.49235505 4 -2.54062227 0.51189162 -0.29552708 1.03157965 0.66279789 0.99295706 5 0.60922674 -0.58659925 0.01509983 -0.18228878 2.09974390 -1.01348568 6 -0.06979795 -1.23783518 0.55391786 0.88632761 -1.08842685 0.28602827 7 1.30693225 0.05841437 -1.16621805 0.42586830 -0.27566626 1.25851116 8 -0.62316734 0.30414603 -0.18446256 0.33951942 -0.13963606 0.05415011 9 0.14145022 1.00802670 0.07073971 -0.86994482 -0.53268266 -0.16444771 10 -0.05753106 -0.31639823 -0.51887891 0.35057350 0.15607134 0.20033499 11 -0.42063131 -0.04186679 0.43208003 -0.41829399 0.07548789 1.26944787 12 > 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/wessaorg/rcomp/tmp/1wt571322607965.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/wessaorg/rcomp/tmp/2nnio1322607965.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/wessaorg/rcomp/tmp/37evd1322607965.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/wessaorg/rcomp/tmp/4lhrm1322607965.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/wessaorg/rcomp/tmp/55u6h1322607965.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/wessaorg/rcomp/createtable file can be downloaded at http://www.wessa.net/cretab > load(file="/var/wessaorg/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/wessaorg/rcomp/tmp/6lj0a1322607965.tab") > > try(system("convert tmp/1wt571322607965.ps tmp/1wt571322607965.png",intern=TRUE)) character(0) > try(system("convert tmp/2nnio1322607965.ps tmp/2nnio1322607965.png",intern=TRUE)) character(0) > try(system("convert tmp/37evd1322607965.ps tmp/37evd1322607965.png",intern=TRUE)) character(0) > try(system("convert tmp/4lhrm1322607965.ps tmp/4lhrm1322607965.png",intern=TRUE)) character(0) > try(system("convert tmp/55u6h1322607965.ps tmp/55u6h1322607965.png",intern=TRUE)) character(0) > > > proc.time() user system elapsed 2.256 0.293 2.564