R version 2.15.2 (2012-10-26) -- "Trick or Treat" Copyright (C) 2012 The R Foundation for Statistical Computing ISBN 3-900051-07-0 Platform: i686-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(54.3,55.9,63.9,64,60.7,67.8,70.5,76.6,76.2,71.8,67.8,69.7,76.7,74.2,75.8,84.3,84.9,84.4,89.4,88.5,76.5,71.4,72.1,75.8,66.6,71.7,75.4,80.9,80.7,85,91.5,87.7,95.3,102.4,114.2,111.7,113.7,118.8,129,136.4,155,166,168.7,145.5,127.3,91.5,69,54,56.3,54.2,59.3,63.4,73.3,86.7,81.3,89.6,85.3,92.4,96.8,93.6,97.6,94.2,99.9,106.4,96,94.9,94.8,95.9,96.2,103.1,106.9,114.2,118.2,123.9,137.1,146.2,136.4,133.2,135.9,127.1,128.5,126.6,132.6,130.9) > par1 = '12' > 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.000000 35.367965 0.000000 6.318842 > m$fitted level slope sea Jan 1 54.30000 0.00000000 0.000000000 Feb 1 55.65290 1.32963888 0.011628435 Mar 1 63.21455 6.38744776 -0.003124094 Apr 1 64.54924 2.24843209 0.005269771 May 1 61.29163 -2.27332987 0.007252018 Jun 1 66.93205 4.22720646 0.007524649 Jul 1 70.55789 3.73322031 0.007496011 Aug 1 76.36680 5.43820904 0.007524306 Sep 1 76.74284 1.28016098 0.007523294 Oct 1 72.40343 -3.33572174 0.007530405 Nov 1 67.91752 -4.28048197 0.007530947 Dec 1 69.09868 0.20568653 0.007530705 Jan 2 75.16771 4.94807802 0.904584866 Feb 2 74.77698 1.09963499 -0.165066804 Mar 2 75.95672 1.16586005 -0.165305268 Apr 2 83.76784 6.61391101 -0.182364740 May 2 85.59877 2.68820556 -0.179358360 Jun 2 84.94302 -0.05859494 -0.179451315 Jul 2 89.11876 3.41972401 -0.179148086 Aug 2 89.05767 0.56053444 -0.179227635 Sep 2 77.94803 -9.02530872 -0.179240957 Oct 2 71.31875 -7.05722297 -0.179245349 Nov 2 71.49304 -1.11729446 -0.179250894 Dec 2 75.42978 3.03406319 -0.179251369 Jan 3 67.03169 -6.26872305 0.799628987 Feb 3 70.79712 1.25071452 0.038797268 Mar 3 75.04506 3.72403571 0.032408059 Apr 3 80.66643 5.28038887 0.028913181 May 3 81.18611 1.37210328 0.031058898 Jun 3 84.73250 3.15809759 0.031102228 Jul 3 91.11785 5.80908240 0.031267912 Aug 3 88.57672 -1.04985547 0.031131103 Sep 3 94.50969 4.68580967 0.031136817 Oct 3 102.05769 7.03679742 0.031133056 Nov 3 113.67128 10.79615257 0.031130540 Dec 3 112.92389 1.31433538 0.031131319 Jan 4 114.07255 1.17916309 -0.354654728 Feb 4 118.44880 3.62591107 0.059910699 Mar 4 128.28100 8.74214832 0.049606118 Apr 4 136.41521 8.24331788 0.050479281 May 4 153.94391 15.86702948 0.047217215 Jun 4 166.33116 13.00870645 0.047163168 Jul 4 169.70118 5.09120796 0.046777506 Aug 4 148.33053 -16.64471670 0.046439610 Sep 4 127.68819 -19.92835441 0.046437060 Oct 4 93.05253 -32.00885594 0.046452125 Nov 4 68.17792 -26.14882259 0.046449068 Dec 4 52.78425 -17.31458593 0.046448503 Jan 5 53.77771 -2.35508049 0.542240140 Feb 5 53.98909 -0.36489518 -0.031104672 Mar 5 58.78110 3.88349738 -0.038116455 Apr 5 63.36335 4.45695298 -0.038938898 May 5 72.79933 8.54526919 -0.040372024 Jun 5 86.21121 12.54273197 -0.040310099 Jul 5 83.04839 -0.35816664 -0.040824920 Aug 5 88.95919 4.79119031 -0.040759339 Sep 5 86.16540 -1.43907774 -0.040763303 Oct 5 91.68430 4.27616429 -0.040769142 Nov 5 96.75445 4.92833983 -0.040769420 Dec 5 94.42936 -1.02959003 -0.040769108 Jan 6 96.64686 1.62611304 0.601632431 Feb 6 94.62831 -1.22584169 -0.076565491 Mar 6 99.34129 3.66449201 -0.083402388 Apr 6 106.14653 6.24239310 -0.086533846 May 6 97.68126 -5.83482457 -0.082948370 Jun 6 94.67535 -3.51113948 -0.082917884 Jul 6 94.51817 -0.75613619 -0.082824773 Aug 6 95.76504 0.88913504 -0.082807027 Sep 6 96.31922 0.61400815 -0.082807175 Oct 6 102.56998 5.24400734 -0.082811181 Nov 6 107.06432 4.62822815 -0.082810958 Dec 6 114.02881 6.54722789 -0.082811043 Jan 7 117.71359 4.20456775 0.796488088 Feb 7 123.78001 5.67109988 -0.062720189 Mar 7 136.41606 11.40450031 -0.069672473 Apr 7 146.42015 10.25494600 -0.068461390 May 7 138.44575 -4.71470197 -0.064607223 Jun 7 133.31035 -5.06026480 -0.064611155 Jul 7 135.20797 0.65504624 -0.064443637 Aug 7 128.01748 -5.78928878 -0.064503919 Sep 7 127.94317 -1.09502920 -0.064501725 Oct 7 126.68251 -1.23107677 -0.064501623 Nov 7 131.95720 4.11272671 -0.064503301 Dec 7 131.46514 0.33037624 -0.064503155 > m$resid Jan Feb Mar Apr May Jun 1 0.00000000 0.24418241 0.86748985 -0.70095063 -0.76054008 1.09306723 2 0.86191692 -0.61666608 0.01096258 0.91121221 -0.66001192 -0.46187089 3 -1.64197256 1.22700706 0.41124988 0.26070161 -0.65710928 0.30031289 4 -0.02354419 0.40260027 0.85281350 -0.08362868 1.28182105 -0.48062398 5 2.58600276 0.32903438 0.70927606 0.09619136 0.68740323 0.67216936 6 0.45684399 -0.47295259 0.81734518 0.43257832 -2.03066485 0.39072541 7 -0.40161703 0.24371956 0.95902271 -0.19295104 -2.51701660 -0.05810607 Jul Aug Sep Oct Nov Dec 1 -0.08306391 0.28669271 -0.69917271 -0.77615727 -0.15886073 0.75434591 2 0.58487400 -0.48077048 -1.61185244 0.33093216 0.99879458 0.69804772 3 0.44575963 -1.15332511 0.96444785 0.39531685 0.63213279 -1.59436055 4 -1.33131801 -3.65487912 -0.55214125 -2.03132739 0.98536027 1.48547036 5 -2.16927188 0.86586046 -1.04761497 0.96101372 0.10966284 -1.00182150 6 0.46325092 0.27665111 -0.04626238 0.77853094 -0.10354282 0.32267839 7 0.96102380 -1.08361007 0.78933628 -0.02287630 0.89855660 -0.63599943 > 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/fisher/rcomp/tmp/1ctlo1354115633.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/fisher/rcomp/tmp/2ambe1354115633.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/fisher/rcomp/tmp/3em531354115633.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/fisher/rcomp/tmp/4hui01354115633.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/fisher/rcomp/tmp/5h2dv1354115633.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/fisher/rcomp/createtable file can be downloaded at http://www.wessa.net/cretab > load(file="/var/fisher/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/fisher/rcomp/tmp/6cllp1354115633.tab") > > try(system("convert tmp/1ctlo1354115633.ps tmp/1ctlo1354115633.png",intern=TRUE)) character(0) > try(system("convert tmp/2ambe1354115633.ps tmp/2ambe1354115633.png",intern=TRUE)) character(0) > try(system("convert tmp/3em531354115633.ps tmp/3em531354115633.png",intern=TRUE)) character(0) > try(system("convert tmp/4hui01354115633.ps tmp/4hui01354115633.png",intern=TRUE)) character(0) > try(system("convert tmp/5h2dv1354115633.ps tmp/5h2dv1354115633.png",intern=TRUE)) character(0) > > > proc.time() user system elapsed 3.440 0.679 4.111