R version 2.12.0 (2010-10-15) Copyright (C) 2010 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(37,30,47,35,30,43,82,40,47,19,52,136,80,42,54,66,81,63,137,72,107,58,36,52,79,77,54,84,48,96,83,66,61,53,30,74,69,59,42,65,70,100,63,105,82,81,75,102,121,98,76,77,63,37,35,23,40,29,37,51,20,28,13,22,25,13,16,13,16,17,9,17,25,14,8,7,10,7,10,3) > 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 66.9018 0.0000 0.0000 308.9231 > m$fitted level slope sea Jan 1 37.000000 0.0000000000 0.0000000000 Feb 1 34.067108 -0.3156477375 -0.3156477363 Mar 1 39.275426 0.1130180120 0.1130180120 Apr 1 37.598825 0.0060909459 0.0060909459 May 1 34.631992 -0.1387018411 -0.1387018411 Jun 1 37.829331 0.0002856516 0.0002856516 Jul 1 54.734562 0.6263475606 0.6263475606 Aug 1 49.275539 0.4201592645 0.4201592645 Sep 1 48.514012 0.3827646184 0.3827646184 Oct 1 37.441837 0.0393364750 0.0393364750 Nov 1 42.949318 0.1961707171 0.1961707171 Dec 1 78.093425 1.1614093467 1.1614093467 Jan 2 82.405129 0.8036404207 -8.8400446320 Feb 2 65.051615 0.0542732866 0.0542732981 Mar 2 60.629143 -0.0882433292 -0.0882433292 Apr 2 62.690874 -0.0332820707 -0.0332820707 May 2 69.659982 0.1173069815 0.1173069815 Jun 2 67.173087 0.0681274596 0.0681274596 Jul 2 93.436458 0.5179943287 0.5179943287 Aug 2 85.531368 0.3829471394 0.3829471394 Sep 2 93.660838 0.5010775977 0.5010775977 Oct 2 80.458127 0.2997695666 0.2997695666 Nov 2 63.928195 0.0595839575 0.0595839575 Dec 2 59.490069 -0.0031260354 -0.0031260354 Jan 3 65.169006 -0.2969582088 3.2665402837 Feb 3 69.932718 -0.1633159500 -0.1633159386 Mar 3 63.696620 -0.2869526044 -0.2869526044 Apr 3 71.361935 -0.1564286797 -0.1564286797 May 3 62.518547 -0.2770175979 -0.2770175979 Jun 3 74.997921 -0.1206463367 -0.1206463367 Jul 3 77.957356 -0.0861220387 -0.0861220387 Aug 3 73.475655 -0.1323633269 -0.1323633269 Sep 3 68.793778 -0.1781178056 -0.1781178056 Oct 3 62.867823 -0.2340777382 -0.2340777382 Nov 3 50.575605 -0.3486754085 -0.3486754085 Dec 3 59.201718 -0.2649593883 -0.2649593883 Jan 4 61.189328 -0.3403712408 3.7440836372 Feb 4 60.250705 -0.3519874709 -0.3519874593 Mar 4 53.165463 -0.4527542008 -0.4527542008 Apr 4 57.526350 -0.3945691722 -0.3945691722 May 4 62.104415 -0.3436056377 -0.3436056377 Jun 4 76.168024 -0.2128682317 -0.2128682317 Jul 4 71.207759 -0.2523643959 -0.2523643959 Aug 4 83.716784 -0.1524647506 -0.1524647506 Sep 4 83.039400 -0.1564029448 -0.1564029448 Oct 4 82.241369 -0.1610743386 -0.1610743386 Nov 4 79.509801 -0.1793862081 -0.1793862081 Dec 4 87.817101 -0.1199088420 -0.1199088420 Jan 5 98.177682 -0.3793096857 4.1724065257 Feb 5 98.025275 -0.3758261367 -0.3758261201 Mar 5 89.551325 -0.4716364840 -0.4716364840 Apr 5 84.714260 -0.5133973314 -0.5133973314 May 5 76.467998 -0.5762115190 -0.5762115190 Jun 5 61.615390 -0.6790513552 -0.6790513552 Jul 5 51.542148 -0.7411892451 -0.7411892451 Aug 5 40.746674 -0.8038613728 -0.8038613728 Sep 5 40.262546 -0.8019484690 -0.8019484690 Oct 5 35.875405 -0.8227936265 -0.8227936265 Nov 5 36.080629 -0.8169367995 -0.8169367995 Dec 5 41.405660 -0.7824643134 -0.7824643134 Jan 6 30.815377 -0.5899092279 6.4890014994 Feb 6 29.592893 -0.5979399281 -0.5979399198 Mar 6 23.164318 -0.6549773977 -0.6549773977 Apr 6 22.563943 -0.6545452069 -0.6545452069 May 6 23.306534 -0.6451458436 -0.6451458436 Jun 6 19.306722 -0.6651804365 -0.6651804365 Jul 6 17.907274 -0.6692111754 -0.6692111754 Aug 6 15.913266 -0.6760714469 -0.6760714469 Sep 6 15.771051 -0.6734154692 -0.6734154692 Oct 6 16.053013 -0.6687918728 -0.6687918728 Nov 6 13.264738 -0.6788522566 -0.6788522566 Dec 6 14.474299 -0.6700135520 -0.6700135520 Jan 7 15.108449 -0.6912306505 7.6035371486 Feb 7 14.522471 -0.6900915928 -0.6900915857 Mar 7 11.892175 -0.7062698244 -0.7062698244 Apr 7 9.884120 -0.7150572269 -0.7150572269 May 7 9.744991 -0.7117506398 -0.7117506398 Jun 7 8.542051 -0.7142558707 -0.7142558707 Jul 7 8.899340 -0.7092287080 -0.7092287080 Aug 7 6.527864 -0.7165904296 -0.7165904296 > m$resid Jan Feb Mar Apr May Jun 1 0.00000000 -0.29155737 0.57810723 -0.19482669 -0.33271936 0.38038993 2 0.53069541 -1.77265794 -0.48610354 0.24498930 0.81704416 -0.30766232 3 0.81616503 0.53927285 -0.68945978 0.92987242 -1.03081304 1.52551249 4 0.30676310 -0.06624776 -0.78005019 0.56963743 0.59475008 1.73293648 5 1.38956318 0.02567529 -0.94906124 -0.52017173 -0.92914822 -1.72302405 6 -1.27958596 -0.07257870 -0.68844343 0.00653550 0.16837916 -0.40577612 7 0.16832630 0.01219335 -0.23028701 -0.15630674 0.06955656 -0.05950699 Jul Aug Sep Oct Nov Dec 1 1.95163537 -0.70835561 -0.13831576 -1.34590443 0.64421930 4.12554663 2 3.11620637 -1.00607626 0.92755315 -1.64344369 -2.02042797 -0.54034739 3 0.36994203 -0.52928193 -0.54865956 -0.69384343 -1.45649857 1.08452209 4 -0.57283476 1.54266392 -0.06352568 -0.07770358 -0.31143553 1.02853430 5 -1.13660581 -1.21823087 0.03877441 -0.43501517 0.12477924 0.74567260 6 -0.08899781 -0.16076480 0.06483022 0.11607004 -0.25757130 0.22953252 7 0.13004200 -0.20193337 > 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/rcomp/tmp/1te4j1322573675.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/rcomp/tmp/2c6vg1322573675.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/rcomp/tmp/35i511322573675.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/rcomp/tmp/483c41322573675.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/rcomp/tmp/595pi1322573675.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/rcomp/createtable file can be downloaded at http://www.wessa.net/cretab > load(file="/var/www/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/rcomp/tmp/62j5w1322573675.tab") > > try(system("convert tmp/1te4j1322573675.ps tmp/1te4j1322573675.png",intern=TRUE)) character(0) > try(system("convert tmp/2c6vg1322573675.ps tmp/2c6vg1322573675.png",intern=TRUE)) character(0) > try(system("convert tmp/35i511322573675.ps tmp/35i511322573675.png",intern=TRUE)) character(0) > try(system("convert tmp/483c41322573675.ps tmp/483c41322573675.png",intern=TRUE)) character(0) > try(system("convert tmp/595pi1322573675.ps tmp/595pi1322573675.png",intern=TRUE)) character(0) > > > proc.time() user system elapsed 2.360 0.220 2.561