R version 2.9.0 (2009-04-17) Copyright (C) 2009 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. 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(16306977,16307888,16307482,16308869,16311019,16312596,16315238,16319511,16327575,16330818,16331930,16334210,16334715,16335459,16334090,16333559,16334600,16336676,16337253,16342333,16348917,16352678,16352972,16357992,16359133,16362938,16365065,16367596,16371278,16374541,16377339,16383275,16393843,16399139,16401009,16405399,16409106,16414307,16418055,16423337,16428686,16434935,16440452,16449092,16464859,16473709,16479291,16485787,16489042,16495231,16501683,16506782,16513615,16520661,16528400,16538542,16554596,16562317,16568499,16574989) > 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.0 2064848.9 202638.1 0.0 > m$fitted level slope sea Jan 1 16306977 0.00000 0.00000 Feb 1 16307770 802.42905 118.35207 Mar 1 16307744 159.00309 -261.92541 Apr 1 16308689 793.07617 179.69731 May 1 16310833 1862.70798 186.07580 Jun 1 16312684 1853.52273 -88.05583 Jul 1 16315143 2334.35982 95.09569 Aug 1 16319238 3732.75074 273.29799 Sep 1 16326957 6899.09117 617.86653 Oct 1 16331475 5008.07321 -657.15749 Nov 1 16332538 1874.89232 -608.15368 Dec 1 16334079 1609.78578 130.75301 Jan 2 16334827 926.09113 -111.74442 Feb 2 16335151 447.84737 307.59993 Mar 2 16334618 -298.04974 -528.11787 Apr 2 16333755 -731.90000 -195.70032 May 2 16334073 71.24427 527.21624 Jun 2 16336231 1662.18099 444.51362 Jul 2 16337524 1380.12465 -271.11666 Aug 2 16342436 4076.61774 -102.61363 Sep 2 16347891 5129.13451 1026.35120 Oct 2 16352699 4883.89714 -20.59683 Nov 2 16354257 2345.22376 -1285.27199 Dec 2 16357468 3005.63222 523.93800 Jan 3 16359406 2191.04633 -272.83284 Feb 3 16362201 2652.35595 737.12978 Mar 3 16365293 2984.19506 -228.32430 Apr 3 16368018 2787.99214 -421.72507 May 3 16371018 2948.79361 259.72688 Jun 3 16373813 2832.51974 728.15542 Jul 3 16378128 3950.47137 -788.79286 Aug 3 16383598 5097.55420 -323.17590 Sep 3 16392197 7740.67530 1645.56965 Oct 3 16398729 6827.78028 410.38820 Nov 3 16403003 4901.53774 -1994.48620 Dec 3 16405095 2782.98676 303.54598 Jan 4 16409192 3774.29150 -86.37478 Feb 4 16413597 4249.49418 710.22369 Mar 4 16418161 4485.10969 -105.54310 Apr 4 16423549 5162.42716 -211.73543 May 4 16428292 4847.38082 394.06123 Jun 4 16434136 5594.66002 799.25895 Jul 4 16441304 6774.24722 -852.04050 Aug 4 16450121 8306.68993 -1029.27228 Sep 4 16462394 11282.40385 2465.03008 Oct 4 16472939 10729.12287 770.07358 Nov 4 16481031 8751.70696 -1740.40896 Dec 4 16486411 6223.40967 -623.72766 Jan 5 16489784 4085.16479 -741.69767 Feb 5 16494443 4515.60734 787.86595 Mar 5 16501481 6399.58274 202.19853 Apr 5 16506931 5690.51396 -149.29486 May 5 16513076 6030.51567 538.82027 Jun 5 16519848 6584.28848 813.45638 Jul 5 16529189 8643.87669 -789.16504 Aug 5 16540342 10518.89452 -1800.17441 Sep 5 16552054 11410.88482 2541.66848 Oct 5 16561520 9956.67740 797.22525 Nov 5 16569739 8658.42531 -1240.08541 Dec 5 16575388 6410.18137 -399.48992 > m$resid Jan Feb Mar Apr May Jun 1 0.000000000 0.596812355 -0.480529304 0.442252698 0.742740018 -0.006397445 2 -0.476802918 -0.337495396 -0.520077636 -0.305053408 0.557250450 1.107259942 3 -0.568275577 0.321476484 0.230791017 -0.137297870 0.111818715 -0.080850396 4 0.691067023 0.330632798 0.163895632 0.472759137 -0.219284536 0.519578141 5 -1.489622481 0.299424907 1.310727204 -0.494327611 0.236716011 0.385112132 Jul Aug Sep Oct Nov Dec 1 0.334609319 0.973170393 2.203506864 -1.315986364 -2.180426185 -0.184491477 2 -0.196343359 1.876454037 0.732463198 -0.170667101 -1.766707401 0.459649120 3 0.778252167 0.798293588 1.839351920 -0.635329031 -1.340477208 -1.474942858 4 0.820976447 1.066542192 2.070802250 -0.385049270 -1.376120142 -1.760401897 5 1.433158192 1.304954962 0.620744531 -1.012012614 -0.903520861 -1.565365304 > 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/rcomp/tmp/1za1m1292851723.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/rcomp/tmp/2ak071292851723.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/rcomp/tmp/3kbia1292851723.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/rcomp/tmp/4kbia1292851723.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/rcomp/tmp/5kbia1292851723.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/rcomp/createtable file can be downloaded at http://www.wessa.net/cretab > load(file="/var/www/html/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/rcomp/tmp/6zlfj1292851723.tab") > > try(system("convert tmp/1za1m1292851723.ps tmp/1za1m1292851723.png",intern=TRUE)) character(0) > try(system("convert tmp/2ak071292851723.ps tmp/2ak071292851723.png",intern=TRUE)) character(0) > try(system("convert tmp/3kbia1292851723.ps tmp/3kbia1292851723.png",intern=TRUE)) character(0) > try(system("convert tmp/4kbia1292851723.ps tmp/4kbia1292851723.png",intern=TRUE)) character(0) > try(system("convert tmp/5kbia1292851723.ps tmp/5kbia1292851723.png",intern=TRUE)) character(0) > > > proc.time() user system elapsed 1.405 0.811 3.192