R version 2.8.0 (2008-10-20) Copyright (C) 2008 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. Natural language support but running in an English locale 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.3156477379 -0.3156477367 Mar 1 39.275426 0.1130180111 0.1130180111 Apr 1 37.598825 0.0060909451 0.0060909451 May 1 34.631992 -0.1387018409 -0.1387018409 Jun 1 37.829331 0.0002856511 0.0002856511 Jul 1 54.734562 0.6263475558 0.6263475558 Aug 1 49.275539 0.4201592597 0.4201592597 Sep 1 48.514012 0.3827646135 0.3827646135 Oct 1 37.441837 0.0393364718 0.0393364718 Nov 1 42.949318 0.1961707141 0.1961707141 Dec 1 78.093425 1.1614093401 1.1614093401 Jan 2 82.405129 0.8036404203 -8.8400446274 Feb 2 65.051615 0.0542732809 0.0542732923 Mar 2 60.629142 -0.0882433346 -0.0882433346 Apr 2 62.690874 -0.0332820758 -0.0332820758 May 2 69.659982 0.1173069759 0.1173069759 Jun 2 67.173087 0.0681274541 0.0681274541 Jul 2 93.436458 0.5179943208 0.5179943208 Aug 2 85.531368 0.3829471312 0.3829471312 Sep 2 93.660838 0.5010775887 0.5010775887 Oct 2 80.458127 0.2997695582 0.2997695582 Nov 2 63.928195 0.0595839504 0.0595839504 Dec 2 59.490069 -0.0031260414 -0.0031260414 Jan 3 65.169006 -0.2969582173 3.2665403767 Feb 3 69.932718 -0.1633159570 -0.1633159455 Mar 3 63.696620 -0.2869526112 -0.2869526112 Apr 3 71.361935 -0.1564286865 -0.1564286865 May 3 62.518547 -0.2770176042 -0.2770176042 Jun 3 74.997922 -0.1206463436 -0.1206463436 Jul 3 77.957356 -0.0861220459 -0.0861220459 Aug 3 73.475655 -0.1323633340 -0.1323633340 Sep 3 68.793778 -0.1781178125 -0.1781178125 Oct 3 62.867823 -0.2340777446 -0.2340777446 Nov 3 50.575604 -0.3486754142 -0.3486754142 Dec 3 59.201718 -0.2649593938 -0.2649593938 Jan 4 61.189328 -0.3403712463 3.7440836979 Feb 4 60.250705 -0.3519874762 -0.3519874646 Mar 4 53.165463 -0.4527542061 -0.4527542061 Apr 4 57.526350 -0.3945691773 -0.3945691773 May 4 62.104415 -0.3436056429 -0.3436056429 Jun 4 76.168024 -0.2128682376 -0.2128682376 Jul 4 71.207759 -0.2523644019 -0.2523644019 Aug 4 83.716784 -0.1524647572 -0.1524647572 Sep 4 83.039400 -0.1564029516 -0.1564029516 Oct 4 82.241369 -0.1610743456 -0.1610743456 Nov 4 79.509801 -0.1793862151 -0.1793862151 Dec 4 87.817101 -0.1199088492 -0.1199088492 Jan 5 98.177682 -0.3793096932 4.1724066086 Feb 5 98.025275 -0.3758261447 -0.3758261281 Mar 5 89.551325 -0.4716364922 -0.4716364922 Apr 5 84.714260 -0.5133973393 -0.5133973393 May 5 76.467998 -0.5762115265 -0.5762115265 Jun 5 61.615390 -0.6790513620 -0.6790513620 Jul 5 51.542148 -0.7411892511 -0.7411892511 Aug 5 40.746674 -0.8038613782 -0.8038613782 Sep 5 40.262546 -0.8019484739 -0.8019484739 Oct 5 35.875405 -0.8227936310 -0.8227936310 Nov 5 36.080629 -0.8169368038 -0.8169368038 Dec 5 41.405660 -0.7824643178 -0.7824643178 Jan 6 30.815377 -0.5899092309 6.4890015317 Feb 6 29.592893 -0.5979399306 -0.5979399223 Mar 6 23.164318 -0.6549774001 -0.6549774001 Apr 6 22.563943 -0.6545452090 -0.6545452090 May 6 23.306534 -0.6451458456 -0.6451458456 Jun 6 19.306722 -0.6651804384 -0.6651804384 Jul 6 17.907274 -0.6692111771 -0.6692111771 Aug 6 15.913266 -0.6760714485 -0.6760714485 Sep 6 15.771051 -0.6734154708 -0.6734154708 Oct 6 16.053013 -0.6687918744 -0.6687918744 Nov 6 13.264738 -0.6788522582 -0.6788522582 Dec 6 14.474299 -0.6700135535 -0.6700135535 Jan 7 15.108449 -0.6912306521 7.6035371659 Feb 7 14.522471 -0.6900915943 -0.6900915873 Mar 7 11.892175 -0.7062698260 -0.7062698260 Apr 7 9.884120 -0.7150572285 -0.7150572285 May 7 9.744991 -0.7117506413 -0.7117506413 Jun 7 8.542051 -0.7142558721 -0.7142558721 Jul 7 8.899340 -0.7092287094 -0.7092287094 Aug 7 6.527864 -0.7165904311 -0.7165904311 > m$resid Jan Feb Mar Apr May Jun 1 0.000000000 -0.291557373 0.578107231 -0.194826690 -0.332719357 0.380389933 2 0.530695403 -1.772657947 -0.486103540 0.244989301 0.817044167 -0.307662326 3 0.816165036 0.539272851 -0.689459782 0.929872419 -1.030813042 1.525512494 4 0.306763103 -0.066247760 -0.780050187 0.569637436 0.594750080 1.732936482 5 1.389563172 0.025675288 -0.949061242 -0.520171726 -0.929148215 -1.723024044 6 -1.279585968 -0.072578700 -0.688443431 0.006535503 0.168379162 -0.405776123 7 0.168326295 0.012193350 -0.230287014 -0.156306743 0.069556557 -0.059506993 Jul Aug Sep Oct Nov Dec 1 1.951635376 -0.708355612 -0.138315764 -1.345904435 0.644219301 4.125546633 2 3.116206371 -1.006076271 0.927553150 -1.643443698 -2.020427969 -0.540347386 3 0.369942030 -0.529281929 -0.548659557 -0.693843429 -1.456498566 1.084522095 4 -0.572834762 1.542663917 -0.063525688 -0.077703579 -0.311435532 1.028534303 5 -1.136605807 -1.218230866 0.038774421 -0.435015164 0.124779247 0.745672601 6 -0.088997804 -0.160764797 0.064830217 0.116070037 -0.257571299 0.229532517 7 0.130042001 -0.201933370 > 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/freestat/rcomp/tmp/1xml11291912287.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/freestat/rcomp/tmp/2qv241291912287.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/freestat/rcomp/tmp/30mjo1291912287.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/freestat/rcomp/tmp/40mjo1291912287.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/freestat/rcomp/tmp/5td0s1291912287.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/freestat/rcomp/createtable file can be downloaded at http://www.wessa.net/cretab > load(file="/var/www/html/freestat/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/freestat/rcomp/tmp/6fwhf1291912287.tab") > > try(system("convert tmp/1xml11291912287.ps tmp/1xml11291912287.png",intern=TRUE)) character(0) > try(system("convert tmp/2qv241291912287.ps tmp/2qv241291912287.png",intern=TRUE)) character(0) > try(system("convert tmp/30mjo1291912287.ps tmp/30mjo1291912287.png",intern=TRUE)) character(0) > try(system("convert tmp/40mjo1291912287.ps tmp/40mjo1291912287.png",intern=TRUE)) character(0) > try(system("convert tmp/5td0s1291912287.ps tmp/5td0s1291912287.png",intern=TRUE)) character(0) > > > proc.time() user system elapsed 2.191 1.211 2.358