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(102.80,118.72,119.01,118.61,120.43,111.83,116.79,131.71,120.57,117.83,130.80,107.46,112.09,129.47,119.72,134.81,135.80,129.27,126.94,153.45,121.86,133.47,135.34,117.10,120.65,132.49,137.60,138.69,125.53,133.09,129.08,145.94,129.07,139.69,142.09,137.29,127.03,137.25,156.87,150.89,139.14,158.30,149.00,158.36,168.06,153.38,173.86,162.47,145.17,168.89,166.64,140.07,128.84,123.40,120.30,129.66,118.12,113.91,131.09,119.14,115.33) > 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 29.88530 0.00000 15.39984 0.00000 > m$fitted level slope sea Jan 1 102.8000 0.0000000 0.00000000 Feb 1 114.7820 0.5447487 3.93796547 Mar 1 118.5741 0.7059074 0.43586319 Apr 1 118.9001 0.6953192 -0.29009132 May 1 119.6135 0.6955644 0.81646126 Jun 1 115.0161 0.6484187 -3.18606411 Jul 1 115.0254 0.6429502 1.76456872 Aug 1 124.4594 0.7228049 7.25062037 Sep 1 123.6151 0.7080741 -3.04505966 Oct 1 119.8907 0.6662152 -2.06072010 Nov 1 125.3283 0.7108975 5.47171745 Dec 1 115.7116 0.6152556 -8.25160679 Jan 2 116.5831 0.6074904 -4.49305842 Feb 2 122.7659 0.6192694 6.70406793 Mar 2 121.9107 0.5882506 -2.19072283 Apr 2 129.0660 0.7272755 5.74403647 May 2 132.2125 0.7626973 3.58746930 Jun 2 133.3468 0.7661747 -4.07677287 Jul 2 131.3258 0.7466721 -4.38584663 Aug 2 138.9195 0.7916739 14.53045740 Sep 2 131.8914 0.7378880 -10.03138660 Oct 2 133.6316 0.7447516 -0.16157126 Nov 2 129.5453 0.7200591 5.79469209 Dec 2 127.0945 0.7192734 -9.99452892 Jan 3 127.1953 0.7229129 -6.54532235 Feb 3 126.7651 0.7208615 5.72494495 Mar 3 135.9955 0.8184175 1.60448211 Apr 3 136.3943 0.8121932 2.29569559 May 3 128.9592 0.7045044 -3.42918158 Jun 3 132.5728 0.7331470 0.51716364 Jul 3 135.0464 0.7461795 -5.96636851 Aug 3 131.8496 0.7211483 14.09037172 Sep 3 135.4763 0.7379560 -6.40629167 Oct 3 137.1084 0.7423187 2.58156141 Nov 3 136.2438 0.7375837 5.84621992 Dec 3 142.2346 0.7393701 -4.94462752 Jan 4 138.9986 0.7418996 -11.96860411 Feb 4 136.7524 0.7332279 0.49755598 Mar 4 146.1444 0.8028693 10.72557976 Apr 4 147.1007 0.8045578 3.78928956 May 4 145.7851 0.7811776 -6.64514786 Jun 4 152.5075 0.8370451 5.79252220 Jul 4 154.1796 0.8433730 -5.17960281 Aug 4 150.0579 0.8127259 8.30213713 Sep 4 163.0938 0.8745020 4.96622362 Oct 4 158.3260 0.8529807 -4.94604237 Nov 4 164.1597 0.8644407 9.70031152 Dec 4 165.9333 0.8653891 -3.46331774 Jan 5 161.6965 0.8593027 -16.52646330 Feb 5 167.2833 0.8754479 1.60674696 Mar 5 162.0340 0.8354337 4.60601130 Apr 5 147.4588 0.7003756 -7.38881829 May 5 140.8215 0.6319448 -11.98146737 Jun 5 127.7102 0.5141984 -4.31019980 Jul 5 124.7227 0.4888277 -4.42266101 Aug 5 125.4445 0.4901982 4.21547458 Sep 5 118.1837 0.4545241 -0.06370772 Oct 5 119.1426 0.4562165 -5.23255268 Nov 5 120.6476 0.4585474 10.44242394 Dec 5 121.0387 0.4584398 -1.89874375 Jan 6 127.4485 0.4704004 -12.11852909 > m$resid Jan Feb Mar Apr May Jun 1 0.000000000 1.501676639 0.495975304 -0.068145177 0.003324129 -0.970778696 2 0.051034284 1.011907981 -0.248991664 1.156236678 0.439502884 0.068017037 3 -0.115157656 -0.209927933 1.498525492 -0.074366429 -1.489915145 0.531001089 4 -0.729692373 -0.542397067 1.547295665 0.027403864 -0.382662049 1.082070861 5 -0.932198632 0.857606664 -1.101522464 -2.767879066 -1.325712407 -2.500337146 6 1.085360452 Jul Aug Sep Oct Nov Dec 1 -0.117094028 1.609676249 -0.286868722 -0.811339238 0.873399486 -1.890502181 2 -0.510570456 1.253875604 -1.431668632 0.183481242 -0.883680351 -0.580831638 3 0.318647599 -0.722314457 0.532171132 0.163656818 -0.293861995 0.962092946 4 0.152748210 -0.909344153 2.238257417 -1.032298420 0.910697186 0.166306047 5 -0.639845405 0.042650525 -1.418571928 0.092236305 0.191717539 -0.012312229 6 > 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/1e76l1259966732.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/2a77t1259966732.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/3gqzi1259966732.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/44sft1259966732.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/5a6qo1259966732.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/6pdve1259966732.tab") > system("convert tmp/1e76l1259966732.ps tmp/1e76l1259966732.png") > system("convert tmp/2a77t1259966732.ps tmp/2a77t1259966732.png") > system("convert tmp/3gqzi1259966732.ps tmp/3gqzi1259966732.png") > system("convert tmp/44sft1259966732.ps tmp/44sft1259966732.png") > system("convert tmp/5a6qo1259966732.ps tmp/5a6qo1259966732.png") > > > proc.time() user system elapsed 1.493 0.800 2.191