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(10.9,10,9.2,9.2,9.5,9.6,9.5,9.1,8.9,9,10.1,10.3,10.2,9.6,9.2,9.3,9.4,9.4,9.2,9,9,9,9.8,10,9.8,9.3,9,9,9.1,9.1,9.1,9.2,8.8,8.3,8.4,8.1,7.7,7.9,7.9,8,7.9,7.6,7.1,6.8,6.5,6.9,8.2,8.7,8.3,7.9,7.5,7.8,8.3,8.4,8.2,7.7,7.2,7.3,8.1,8.5) > 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.0000000000 0.1607481629 0.0008262046 0.0000000000 > m$fitted level slope sea Jan 1 10.900000 0.000000000 0.000000e+00 Feb 1 10.006864 -0.893703334 -6.863933e-03 Mar 1 9.193475 -0.815149386 6.525104e-03 Apr 1 9.185281 -0.030578051 1.471936e-02 May 1 9.497454 0.302822348 2.546057e-03 Jun 1 9.605338 0.113190193 -5.338419e-03 Jul 1 9.503155 -0.096318749 -3.155085e-03 Aug 1 9.104686 -0.390241136 -4.686163e-03 Sep 1 8.895208 -0.214400853 4.792497e-03 Oct 1 8.995008 0.091244988 4.991748e-03 Nov 1 10.082837 1.060691415 1.716318e-02 Dec 1 10.320254 0.259834811 -2.025372e-02 Jan 2 10.203188 -0.106739137 -3.187889e-03 Feb 2 9.606446 -0.583545988 -6.446039e-03 Mar 2 9.199129 -0.415052931 8.711569e-04 Apr 2 9.282231 0.060182915 1.776937e-02 May 2 9.396844 0.112109689 3.156391e-03 Jun 2 9.408928 0.016671325 -8.927668e-03 Jul 2 9.199966 -0.198611874 3.404688e-05 Aug 2 9.004655 -0.195462147 -4.655257e-03 Sep 2 8.978978 -0.033466092 2.102163e-02 Oct 2 9.032773 0.049791665 -3.277326e-02 Nov 2 9.745057 0.681893668 5.494289e-02 Dec 2 10.036128 0.309002380 -3.612756e-02 Jan 3 9.803109 -0.208058299 -3.108662e-03 Feb 3 9.307446 -0.482573050 -7.446320e-03 Mar 3 9.009526 -0.308339324 -9.526329e-03 Apr 3 8.976922 -0.048147400 2.307777e-02 May 3 9.094113 0.107797804 5.886545e-03 Jun 3 9.104205 0.015625212 -4.205341e-03 Jul 3 9.097779 -0.005177213 2.220567e-03 Aug 3 9.205782 0.101591825 -5.781741e-03 Sep 3 8.784378 -0.391780570 1.562221e-02 Oct 3 8.369136 -0.413913168 -6.913579e-02 Nov 3 8.325530 -0.064577075 7.447013e-02 Dec 3 8.126998 -0.190941601 -2.699837e-02 Jan 4 7.695290 -0.418029543 4.710068e-03 Feb 4 7.883674 0.154181344 1.632610e-02 Mar 4 7.928428 0.051661599 -2.842782e-02 Apr 4 7.990096 0.061044694 9.904375e-03 May 4 7.893690 -0.086495313 6.309921e-03 Jun 4 7.607373 -0.273778118 -7.373376e-03 Jul 4 7.127907 -0.466564351 -2.790678e-02 Aug 4 6.769976 -0.364746080 3.002360e-02 Sep 4 6.468983 -0.304993457 3.101686e-02 Oct 4 6.963943 0.444770777 -6.394334e-02 Nov 4 8.100087 1.092771501 9.991328e-02 Dec 4 8.720299 0.649881733 -2.029946e-02 Jan 5 8.373373 -0.284241898 -7.337275e-02 Feb 5 7.880350 -0.479931048 1.965010e-02 Mar 5 7.524156 -0.364474447 -2.415621e-02 Apr 5 7.760586 0.196821075 3.941389e-02 May 5 8.276755 0.494896693 2.324533e-02 Jun 5 8.412750 0.159842838 -1.275004e-02 Jul 5 8.246627 -0.144476337 -4.662740e-02 Aug 5 7.661113 -0.556222340 3.888715e-02 Sep 5 7.189124 -0.477583864 1.087589e-02 Oct 5 7.393795 0.159356857 -9.379471e-02 Nov 5 7.986242 0.563684299 1.137579e-01 Dec 5 8.472987 0.491860585 2.701250e-02 > m$resid Jan Feb Mar Apr May Jun 1 0.000000000 -2.237246511 0.198063229 1.956005606 0.831574902 -0.472975904 2 -0.914533314 -1.192633213 0.423495729 1.184434295 0.129523187 -0.238039611 3 -1.290020970 -0.685783396 0.436606586 0.648553124 0.388980473 -0.229894430 4 -0.566605393 1.428172029 -0.256407404 0.023394412 -0.367995665 -0.467119686 5 -2.330819063 -0.488182429 0.288448811 1.399766343 0.743416939 -0.835698086 Jul Aug Sep Oct Nov Dec 1 -0.522552051 -0.733093986 0.438576505 0.762334339 2.417969424 -1.997476838 2 -0.536954217 0.007855973 0.404046576 0.207659448 1.576573427 -0.930057139 3 -0.051884925 0.266300706 -1.230557310 -0.055202583 0.871306220 -0.315175885 4 -0.480842524 0.253952669 0.149033529 1.870043551 1.616231020 -1.104648967 5 -0.759024030 -1.026967121 0.196138158 1.588642189 1.008466422 -0.179142886 > 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/1cmvq1259934900.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/2vc021259934900.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/3hixw1259934900.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/4r0iq1259934900.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/5jw1u1259934900.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/6gq8u1259934900.tab") > system("convert tmp/1cmvq1259934900.ps tmp/1cmvq1259934900.png") > system("convert tmp/2vc021259934900.ps tmp/2vc021259934900.png") > system("convert tmp/3hixw1259934900.ps tmp/3hixw1259934900.png") > system("convert tmp/4r0iq1259934900.ps tmp/4r0iq1259934900.png") > system("convert tmp/5jw1u1259934900.ps tmp/5jw1u1259934900.png") > > > proc.time() user system elapsed 1.362 0.825 1.801