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(108.8,128.4,121.1,119.5,128.7,108.7,105.5,119.8,111.3,110.6,120.1,97.5,107.7,127.3,117.2,119.8,116.2,111,112.4,130.6,109.1,118.8,123.9,101.6,112.8,128,129.6,125.8,119.5,115.7,113.6,129.7,112,116.8,127,112.1,114.2,121.1,131.6,125,120.4,117.7,117.5,120.6,127.5,112.3,124.5,115.2,104.7,130.9,129.2,113.5,125.6,107.6,107,121.6,110.7,106.3,118.6,104.6) > 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 7.6905088 0.0000000 16.4615192 0.0617548 > m$fitted level slope sea Jan 1 108.8000 0.0000000000 0.0000000 Feb 1 117.2112 0.0955901789 11.1621325 Mar 1 122.6133 0.4467422859 -1.5225623 Apr 1 122.4317 0.4066493826 -2.9303194 May 1 124.7643 0.4970153529 3.9305757 Jun 1 119.4284 0.3144275309 -10.7114641 Jul 1 112.6531 0.1584652668 -7.1321140 Aug 1 113.2725 0.1667286266 6.5261147 Sep 1 112.7148 0.1545722233 -1.4126261 Oct 1 111.7485 0.1357480879 -1.1451752 Nov 1 114.5740 0.1815367604 5.5180290 Dec 1 109.0289 0.0835286661 -11.5119123 Jan 2 109.2404 0.0814538546 -1.5408322 Feb 2 112.8831 0.0563866848 14.4060152 Mar 2 116.0725 0.1022301661 1.1188939 Apr 2 119.3240 0.1871965071 0.4678771 May 2 115.6740 0.0745110154 0.5360297 Jun 2 115.5747 0.0699802926 -4.5742028 Jul 2 117.1098 0.1010200408 -4.7139271 Aug 2 119.8384 0.1458518486 10.7541341 Sep 2 117.0491 0.1043995642 -7.9407548 Oct 2 117.3847 0.1071051860 1.4146226 Nov 2 116.7021 0.1001137548 7.2001834 Dec 2 115.7702 0.0948073598 -14.1672090 Jan 3 116.2641 0.0955408814 -3.4652830 Feb 3 116.3537 0.0955183354 11.6463125 Mar 3 121.1579 0.1419198929 8.4289502 Apr 3 123.1795 0.1716936226 2.6153458 May 3 122.4185 0.1538511674 -2.9160181 Jun 3 122.0447 0.1435870387 -6.3432350 Jul 3 120.9568 0.1217355211 -7.3533733 Aug 3 119.3200 0.0951266422 10.3849044 Sep 3 119.0605 0.0907519078 -7.0595167 Oct 3 117.4033 0.0738979531 -0.5983699 Nov 3 117.7054 0.0755404095 9.2939526 Dec 3 121.1836 0.0937628455 -9.0933012 Jan 4 121.4284 0.0944847539 -7.2287976 Feb 4 118.4417 0.0758275817 2.6670281 Mar 4 119.7220 0.0864262258 11.8746513 Apr 4 120.7093 0.0970814467 4.2881796 May 4 121.5779 0.1078836883 -1.1800132 Jun 4 122.1826 0.1152284657 -4.4839702 Jul 4 122.8063 0.1224632561 -5.3077639 Aug 4 118.2628 0.0630692139 2.3502762 Sep 4 122.7828 0.1110443003 4.7046099 Oct 4 120.8354 0.0930351472 -8.5296085 Nov 4 119.0201 0.0796031427 5.4852988 Dec 4 120.2534 0.0864803833 -5.0567076 Jan 5 117.1111 0.0678060931 -12.4019560 Feb 5 120.7864 0.0915284173 10.1034808 Mar 5 120.8159 0.0910275097 8.3842820 Apr 5 116.8368 0.0512256512 -3.3254406 May 5 119.4996 0.0802766178 6.0931808 Jun 5 117.2935 0.0534024354 -9.6871527 Jul 5 114.3089 0.0182361599 -7.3004597 Aug 5 116.1016 0.0372622941 5.4933919 Sep 5 112.3092 0.0010210155 -1.5984115 Oct 5 111.8956 -0.0023424688 -5.5944184 Nov 5 112.3193 0.0006253244 6.2794901 Dec 5 111.0985 -0.0070249789 -6.4950473 > m$resid Jan Feb Mar Apr May Jun 1 0.000000000 2.910424152 1.275843874 -0.185107233 0.652755826 -2.070628326 2 0.048147476 1.330820135 1.069177125 1.036830124 -1.294380452 -0.060591677 3 0.144212801 -0.002139232 1.650053921 0.647216238 -0.321529572 -0.184297305 4 0.054135274 -1.098002313 0.424949908 0.315284817 0.269681697 0.174596772 5 -1.154485809 1.284742151 -0.021959638 -1.435786308 0.920351498 -0.807659119 Jul Aug Sep Oct Nov Dec 1 -2.544767190 0.165735950 -0.260287108 -0.402205333 0.964273439 -2.052057928 2 0.520614000 0.941596110 -1.054871037 0.083107448 -0.283758825 -0.371235483 3 -0.435876160 -0.627821540 -0.127128210 -0.627323265 0.081901535 1.221551011 4 0.180059858 -1.662411812 1.593721335 -0.737249491 -0.683842588 0.413279777 5 -1.077608921 0.631925100 -1.367634153 -0.148292488 0.152464753 -0.436941680 > 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/1nrsj1260522912.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/2y4xo1260522912.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/32vqk1260522912.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/4epwk1260522912.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/550wt1260522912.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/6uf531260522912.tab") > system("convert tmp/1nrsj1260522912.ps tmp/1nrsj1260522912.png") > system("convert tmp/2y4xo1260522912.ps tmp/2y4xo1260522912.png") > system("convert tmp/32vqk1260522912.ps tmp/32vqk1260522912.png") > system("convert tmp/4epwk1260522912.ps tmp/4epwk1260522912.png") > system("convert tmp/550wt1260522912.ps tmp/550wt1260522912.png") > > > proc.time() user system elapsed 1.727 0.835 6.421