R version 2.12.0 (2010-10-15) Copyright (C) 2010 The R Foundation for Statistical Computing ISBN 3-900051-07-0 Platform: x86_64-redhat-linux-gnu (64-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(31245,30951,30872,30752,30967,30781,30681,31356,31434,31594,31949,32396,32441,32447,32288,32418,32346,32091,31855,31683,31615,31840,31536,31383,31638,31626,31720,31472,31372,31419,31341,31171,31036,30532,30666,30571,30173,30032,29874,30018,29911,29963,30050,29901,29544,29451,29293,29334,29389,29563) > 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 44643.8325 126.8165 0.0000 0.0000 > m$fitted level slope sea Jan 1 31245.00 0.0000000 0.000000 Feb 1 30966.27 -15.5442802 -15.273942 Mar 1 30887.52 -16.1346587 -15.516608 Apr 1 30767.91 -17.3809076 -15.909043 May 1 30982.04 -13.9829550 -15.043877 Jun 1 30796.67 -16.9353258 -15.673313 Jul 1 30696.97 -18.5643541 -15.971298 Aug 1 31369.54 -3.3342627 -13.537849 Sep 1 31447.26 -1.3660026 -13.259385 Oct 1 31606.72 2.8824969 -12.721461 Nov 1 31960.58 12.8623014 -11.580922 Dec 1 32406.22 25.9905887 -10.217237 Jan 2 32361.24 26.3057226 79.756778 Feb 2 32452.70 29.5342131 -5.697055 Mar 2 32293.94 23.1125478 -5.939577 Apr 2 32423.81 26.9164977 -5.806976 May 2 32351.93 23.2559089 -5.925147 Jun 2 32097.24 12.5930009 -6.244830 Jul 2 31861.52 2.7651558 -6.519143 Aug 2 31689.70 -4.3393798 -6.704149 Sep 2 31621.77 -6.9927397 -6.768732 Oct 2 31846.54 2.8954430 -6.543411 Nov 2 31542.83 -10.4503340 -6.828518 Dec 2 31389.96 -16.7618423 -6.955084 Jan 3 31527.83 -11.4293163 110.165467 Feb 3 31634.24 -5.2802657 -8.243689 Mar 3 31728.20 -0.7026506 -8.197021 Apr 3 31480.31 -12.2426510 -8.307837 May 3 31380.35 -16.3821605 -8.345307 Jun 3 31427.32 -13.3633545 -8.319533 Jul 3 31349.34 -16.4687895 -8.344554 Aug 3 31179.40 -23.9029554 -8.401108 Sep 3 31044.44 -29.3202906 -8.440036 Oct 3 30540.60 -52.6131253 -8.598198 Nov 3 30674.54 -43.4038667 -8.539088 Dec 3 30579.55 -45.9631069 -8.554620 Jan 4 30136.43 -64.2097513 36.566342 Feb 4 30035.72 -66.1318609 -3.723530 Mar 4 29877.74 -70.7431145 -3.738683 Apr 4 30021.71 -59.9278245 -3.705047 May 4 29914.71 -62.3057582 -3.712048 Jun 4 29966.70 -56.5156469 -3.695909 Jul 4 30053.68 -49.2280955 -3.676676 Aug 4 29904.69 -54.3055567 -3.689366 Sep 4 29547.73 -69.7403941 -3.725905 Oct 4 29454.73 -70.9285469 -3.728569 Nov 4 29296.74 -75.3834673 -3.738032 Dec 4 29337.73 -69.4202664 -3.726031 Jan 5 29340.18 -65.8489638 48.824269 Feb 5 29564.94 -50.6158566 -1.937521 > m$resid Jan Feb Mar Apr May Jun Jul 1 0.0000000 -0.7916612 -0.2989237 -0.4886171 1.0917446 -0.8071076 -0.3892555 2 -0.3981492 0.2588059 -0.8769691 0.4967950 -0.4594091 -1.2914537 -1.1530528 3 0.7854012 0.4970256 0.4589145 -1.1427773 -0.4054250 0.2927450 -0.2984739 4 -1.9508707 -0.1579544 -0.4237421 0.9904148 -0.2170837 0.5271079 0.6617626 5 0.3473751 1.2756880 Aug Sep Oct Nov Dec 1 3.2462691 0.3802536 0.7535942 1.6427200 2.0233868 2 -0.8101440 -0.2949480 1.0743238 -1.4205459 -0.6595561 3 -0.7088238 -0.5128216 -2.1907563 0.8611363 -0.2380616 4 -0.4600295 -1.3955930 -0.1072348 -0.4014123 0.5365242 5 > 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/wessaorg/rcomp/tmp/1d08k1293639774.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/wessaorg/rcomp/tmp/2d08k1293639774.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/wessaorg/rcomp/tmp/3o9q51293639774.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/wessaorg/rcomp/tmp/4z0p81293639774.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/wessaorg/rcomp/tmp/5z0p81293639774.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/wessaorg/rcomp/createtable file can be downloaded at http://www.wessa.net/cretab > load(file="/var/www/wessaorg/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/wessaorg/rcomp/tmp/6nkqt1293639775.tab") > > try(system("convert tmp/1d08k1293639774.ps tmp/1d08k1293639774.png",intern=TRUE)) character(0) > try(system("convert tmp/2d08k1293639774.ps tmp/2d08k1293639774.png",intern=TRUE)) character(0) > try(system("convert tmp/3o9q51293639774.ps tmp/3o9q51293639774.png",intern=TRUE)) character(0) > try(system("convert tmp/4z0p81293639774.ps tmp/4z0p81293639774.png",intern=TRUE)) character(0) > try(system("convert tmp/5z0p81293639774.ps tmp/5z0p81293639774.png",intern=TRUE)) character(0) > > > proc.time() user system elapsed 1.65 0.16 1.93