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(-5,-1,-2,-5,-4,-6,-2,-2,-2,-2,2,1,-8,-1,1,-1,2,2,1,-1,-2,-2,-1,-8,-4,-6,-3,-3,-7,-9,-11,-13,-11,-9,-17,-22,-25,-20,-24,-24,-22,-19,-18,-17,-11,-11,-12,-10,-15,-15,-15,-13,-8,-13,-9,-7,-4,-4,-2,0) > 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 4.48673122 0.04721141 0.00000000 2.40674812 > m$fitted level slope sea Jan 1 -5.0000000 0.00000000 0.0000000 Feb 1 -1.7738294 0.21127491 0.2000189 Mar 1 -2.0247942 0.19354403 0.1890231 Apr 1 -4.2606995 0.10400284 0.1599449 May 1 -4.1590561 0.10390258 0.1599250 Jun 1 -5.5896436 0.02803981 0.1489405 Jul 1 -3.0655919 0.16809565 0.1651745 Aug 1 -2.3603530 0.20145149 0.1684045 Sep 1 -2.1658893 0.20098104 0.1683656 Oct 1 -2.1148485 0.19021284 0.1675944 Nov 1 0.8495366 0.40004868 0.1807273 Dec 1 0.9313090 0.37494785 0.1793450 Jan 2 -2.3430512 0.20566284 -4.4167470 Feb 2 -1.5641538 0.25936732 0.4173219 Mar 2 0.0982212 0.38201053 0.4296079 Apr 2 -0.9395305 0.25845392 0.4250593 May 2 0.9985304 0.40610780 0.4277184 Jun 2 1.5295586 0.41721844 0.4278516 Jul 2 0.9218811 0.32511266 0.4269666 Aug 2 -0.7475065 0.14431313 0.4254642 Sep 2 -1.9628843 0.02018193 0.4245481 Oct 2 -2.3023338 -0.01284025 0.4243297 Nov 2 -1.6501254 0.04851347 0.4246944 Dec 2 -6.6963022 -0.42329191 0.4221703 Jan 3 -2.4724387 -0.02040723 -3.0598015 Feb 3 -5.3986339 -0.29630562 0.2270768 Mar 3 -3.8535569 -0.12381609 0.2362322 Apr 3 -3.4241159 -0.07207189 0.2369822 May 3 -6.2909567 -0.33360848 0.2360085 Jun 3 -8.5761869 -0.51641307 0.2358723 Jul 3 -10.6945181 -0.66659450 0.2358810 Aug 3 -12.7624135 -0.79806728 0.2359146 Sep 3 -11.8229306 -0.63494482 0.2358676 Oct 3 -10.0494680 -0.40872578 0.2358052 Nov 3 -15.5244364 -0.88479158 0.2359258 Dec 3 -20.7647121 -1.29420838 0.2360199 Jan 4 -22.2796429 -1.31455686 -2.6473856 Feb 4 -21.0382648 -1.07389376 0.2739941 Mar 4 -23.7262199 -1.22569160 0.2682877 Apr 4 -24.4412564 -1.17764973 0.2687175 May 4 -23.1148409 -0.94203583 0.2689755 Jun 4 -20.4775030 -0.60521668 0.2686515 Jul 4 -18.9788587 -0.40723078 0.2683717 Aug 4 -17.8028295 -0.25822834 0.2681589 Sep 4 -12.9825513 0.21973857 0.2675226 Oct 4 -11.6449152 0.32495448 0.2673950 Nov 4 -12.0282674 0.25828675 0.2674682 Dec 4 -10.6467052 0.36401552 0.2673634 Jan 5 -11.4609627 0.25413497 -3.1478321 Feb 5 -14.3034362 -0.03669329 0.2539278 Mar 5 -15.0225581 -0.10087381 0.2520282 Apr 5 -13.7252567 0.03072322 0.2529128 May 5 -9.6269456 0.41363395 0.2530984 Jun 5 -12.2334992 0.12931611 0.2534177 Jul 5 -9.9728507 0.32995683 0.2531268 Aug 5 -7.8562232 0.49815041 0.2528862 Sep 5 -5.0365543 0.71669247 0.2525967 Oct 5 -4.2695733 0.72142650 0.2525910 Nov 5 -2.5795778 0.81260491 0.2524917 Dec 5 -0.6347450 0.91918976 0.2523868 > m$resid Jan Feb Mar Apr May Jun 1 0.000000000 0.976560021 -0.208563075 -1.120081036 -0.001082880 -0.699372217 2 -1.943951226 0.211882082 0.608518389 -0.620712732 0.733360810 0.054466664 3 2.223186822 -1.141476907 0.795613255 0.240055167 -1.211977080 -0.846044928 4 -0.102254673 1.031951444 -0.697786862 0.221379194 1.085100431 1.550756958 5 -0.537337055 -1.269269134 -0.295170460 0.606034797 1.762454979 -1.308418365 Jul Aug Sep Oct Nov Dec 1 1.129934682 0.241672950 -0.003126877 -0.066775258 1.230366380 -0.140656663 2 -0.446370303 -0.867887544 -0.591226694 -0.156282871 0.288863270 -2.211967718 3 -0.694332627 -0.607321959 0.753005199 1.043686788 -2.195382142 -1.887329360 4 0.911444138 0.685901720 2.200127991 0.484303383 -0.306860201 0.486641430 5 0.923311794 0.774002599 1.005714738 0.021786000 0.419608067 0.490514508 > 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/1qynp1291727317.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/html/rcomp/tmp/21p5a1291727317.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/html/rcomp/tmp/31p5a1291727317.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/html/rcomp/tmp/4uymv1291727317.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/www/html/rcomp/tmp/5uymv1291727317.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/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/6iimh1291727318.tab") > > try(system("convert tmp/1qynp1291727317.ps tmp/1qynp1291727317.png",intern=TRUE)) character(0) > try(system("convert tmp/21p5a1291727317.ps tmp/21p5a1291727317.png",intern=TRUE)) character(0) > try(system("convert tmp/31p5a1291727317.ps tmp/31p5a1291727317.png",intern=TRUE)) character(0) > try(system("convert tmp/4uymv1291727317.ps tmp/4uymv1291727317.png",intern=TRUE)) character(0) > try(system("convert tmp/5uymv1291727317.ps tmp/5uymv1291727317.png",intern=TRUE)) character(0) > > > proc.time() user system elapsed 1.387 0.838 3.995