R version 2.12.1 (2010-12-16) 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.3156477376 -0.3156477364 Mar 1 39.275426 0.1130180118 0.1130180118 Apr 1 37.598825 0.0060909457 0.0060909457 May 1 34.631992 -0.1387018411 -0.1387018411 Jun 1 37.829331 0.0002856515 0.0002856515 Jul 1 54.734562 0.6263475595 0.6263475595 Aug 1 49.275539 0.4201592633 0.4201592633 Sep 1 48.514012 0.3827646172 0.3827646172 Oct 1 37.441837 0.0393364742 0.0393364742 Nov 1 42.949318 0.1961707164 0.1961707164 Dec 1 78.093425 1.1614093451 1.1614093451 Jan 2 82.405129 0.8036404207 -8.8400446315 Feb 2 65.051615 0.0542732851 0.0542732966 Mar 2 60.629143 -0.0882433306 -0.0882433306 Apr 2 62.690874 -0.0332820721 -0.0332820721 May 2 69.659982 0.1173069801 0.1173069801 Jun 2 67.173087 0.0681274581 0.0681274581 Jul 2 93.436458 0.5179943267 0.5179943267 Aug 2 85.531368 0.3829471373 0.3829471373 Sep 2 93.660838 0.5010775954 0.5010775954 Oct 2 80.458127 0.2997695645 0.2997695645 Nov 2 63.928195 0.0595839557 0.0595839557 Dec 2 59.490069 -0.0031260370 -0.0031260370 Jan 3 65.169006 -0.2969582110 3.2665403080 Feb 3 69.932718 -0.1633159518 -0.1633159404 Mar 3 63.696620 -0.2869526062 -0.2869526062 Apr 3 71.361935 -0.1564286815 -0.1564286815 May 3 62.518547 -0.2770175995 -0.2770175995 Jun 3 74.997921 -0.1206463385 -0.1206463385 Jul 3 77.957356 -0.0861220405 -0.0861220405 Aug 3 73.475655 -0.1323633288 -0.1323633288 Sep 3 68.793778 -0.1781178074 -0.1781178074 Oct 3 62.867823 -0.2340777398 -0.2340777398 Nov 3 50.575605 -0.3486754100 -0.3486754100 Dec 3 59.201718 -0.2649593898 -0.2649593898 Jan 4 61.189328 -0.3403712422 3.7440836531 Feb 4 60.250705 -0.3519874723 -0.3519874607 Mar 4 53.165463 -0.4527542022 -0.4527542022 Apr 4 57.526350 -0.3945691735 -0.3945691735 May 4 62.104415 -0.3436056391 -0.3436056391 Jun 4 76.168024 -0.2128682332 -0.2128682332 Jul 4 71.207759 -0.2523643975 -0.2523643975 Aug 4 83.716784 -0.1524647523 -0.1524647523 Sep 4 83.039400 -0.1564029465 -0.1564029465 Oct 4 82.241369 -0.1610743404 -0.1610743404 Nov 4 79.509801 -0.1793862099 -0.1793862099 Dec 4 87.817101 -0.1199088439 -0.1199088439 Jan 5 98.177682 -0.3793096876 4.1724065473 Feb 5 98.025275 -0.3758261388 -0.3758261222 Mar 5 89.551325 -0.4716364861 -0.4716364861 Apr 5 84.714260 -0.5133973334 -0.5133973334 May 5 76.467998 -0.5762115210 -0.5762115210 Jun 5 61.615390 -0.6790513570 -0.6790513570 Jul 5 51.542148 -0.7411892467 -0.7411892467 Aug 5 40.746674 -0.8038613742 -0.8038613742 Sep 5 40.262546 -0.8019484703 -0.8019484703 Oct 5 35.875405 -0.8227936277 -0.8227936277 Nov 5 36.080629 -0.8169368006 -0.8169368006 Dec 5 41.405660 -0.7824643146 -0.7824643146 Jan 6 30.815377 -0.5899092287 6.4890015080 Feb 6 29.592893 -0.5979399287 -0.5979399205 Mar 6 23.164318 -0.6549773983 -0.6549773983 Apr 6 22.563943 -0.6545452074 -0.6545452074 May 6 23.306534 -0.6451458441 -0.6451458441 Jun 6 19.306722 -0.6651804370 -0.6651804370 Jul 6 17.907274 -0.6692111759 -0.6692111759 Aug 6 15.913266 -0.6760714473 -0.6760714473 Sep 6 15.771051 -0.6734154697 -0.6734154697 Oct 6 16.053013 -0.6687918732 -0.6687918732 Nov 6 13.264738 -0.6788522570 -0.6788522570 Dec 6 14.474299 -0.6700135524 -0.6700135524 Jan 7 15.108449 -0.6912306509 7.6035371533 Feb 7 14.522471 -0.6900915932 -0.6900915861 Mar 7 11.892175 -0.7062698248 -0.7062698248 Apr 7 9.884120 -0.7150572273 -0.7150572273 May 7 9.744991 -0.7117506402 -0.7117506402 Jun 7 8.542051 -0.7142558711 -0.7142558711 Jul 7 8.899340 -0.7092287083 -0.7092287083 Aug 7 6.527864 -0.7165904300 -0.7165904300 > 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.81704417 -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.27958597 -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/1wnoy1322744561.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/2nd8v1322744561.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/39i501322744561.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/4s1x91322744561.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/5784b1322744561.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/6ykfl1322744561.tab") > > try(system("convert tmp/1wnoy1322744561.ps tmp/1wnoy1322744561.png",intern=TRUE)) character(0) > try(system("convert tmp/2nd8v1322744561.ps tmp/2nd8v1322744561.png",intern=TRUE)) character(0) > try(system("convert tmp/39i501322744561.ps tmp/39i501322744561.png",intern=TRUE)) character(0) > try(system("convert tmp/4s1x91322744561.ps tmp/4s1x91322744561.png",intern=TRUE)) character(0) > try(system("convert tmp/5784b1322744561.ps tmp/5784b1322744561.png",intern=TRUE)) character(0) > > > proc.time() user system elapsed 2.696 0.316 2.986