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(10570,10297,10635,10872,10296,10383,10431,10574,10653,10805,10872,10625,10407,10463,10556,10646,10702,11353,11346,11451,11964,12574,13031,13812,14544,14931,14886,16005,17064,15168,16050,15839,15137,14954,15648,15305,15579,16348,15928,16171,15937,15713,15594,15683,16438,17032,17696,17745,19394,20148,20108,18584,18441,18391,19178,18079,18483,19644,19195,19650,20830,23595,22937,21814,21928,21777,21383,21467,22052,22680,24320,24977,25204) > 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 467254.0204 871.5329 903.4623 0.0000 > m$fitted level slope sea Jan 1 10570.00 0.00000000 0.0000000 Feb 1 10311.43 -14.33001689 -14.4309339 Mar 1 10645.43 -11.55036410 -10.4340168 Apr 1 10882.93 -9.24858912 -10.9297069 May 1 10313.53 -15.45247856 -17.5339427 Jun 1 10394.86 -14.21530915 -11.8579236 Jul 1 10444.06 -13.29923346 -13.0648803 Aug 1 10585.91 -10.80859595 -11.9057150 Sep 1 10665.13 -9.22316569 -12.1256982 Oct 1 10816.03 -6.16310036 -11.0338840 Nov 1 10883.40 -4.65198392 -11.4040079 Dec 1 10638.84 -9.91323577 -13.8389634 Jan 2 10383.94 -5.96314347 23.0555156 Feb 2 10464.55 -2.36986940 -1.5467515 Mar 2 10555.00 -0.07204888 1.0046164 Apr 2 10640.95 2.11063090 5.0518406 May 2 10705.37 3.76786453 -3.3663934 Jun 2 11343.40 21.32496664 9.6026926 Jul 2 11347.81 20.83971061 -1.8142531 Aug 2 11448.19 23.19880877 2.8079960 Sep 2 11957.53 38.05687072 6.4665090 Oct 2 12566.95 56.00435464 7.0464905 Nov 2 13024.38 68.94260342 6.6244681 Dec 2 13802.36 92.21532031 9.6355171 Jan 3 14461.06 104.00092395 82.9431154 Feb 3 14933.02 119.56009198 -2.0239421 Mar 3 14892.23 113.93097431 -6.2272792 Apr 3 15983.77 148.32056130 21.2301912 May 3 17065.87 181.72868108 -1.8742391 Jun 3 15199.03 107.40077434 -31.0330664 Jul 3 16031.84 134.05640944 18.1556002 Aug 3 15851.32 122.36180069 -12.3175930 Sep 3 15151.32 91.46344895 -14.3172989 Oct 3 14954.32 80.51865768 -0.3242140 Nov 3 15641.16 103.74742326 6.8385562 Dec 3 15316.57 87.28141837 -11.5652746 Jan 4 15532.78 91.61101708 46.2194085 Feb 4 16328.13 121.77757098 19.8688530 Mar 4 15954.50 102.12300081 -26.4988769 Apr 4 16164.80 106.40020902 6.1966574 May 4 15906.83 91.90203917 30.1708987 Jun 4 15755.08 82.15734597 -42.0777473 Jul 4 15574.18 71.58904912 19.8156488 Aug 4 15676.39 72.82432431 6.6075075 Sep 4 16425.92 100.22926099 12.0750474 Oct 4 17032.48 120.80715767 -0.4770064 Nov 4 17668.36 141.81891065 27.6406366 Dec 4 17757.42 139.66897019 -12.4228853 Jan 5 19274.87 192.86043811 119.1311149 Feb 5 20115.83 220.62397475 32.1749030 Mar 5 20157.24 213.20511002 -49.2353552 Apr 5 18621.14 141.06650173 -37.1449944 May 5 18410.78 126.53932694 30.2209288 Jun 5 18433.85 122.25405426 -42.8517844 Jul 5 19146.83 146.75996799 31.1690879 Aug 5 18122.96 98.12460669 -43.9570800 Sep 5 18471.69 108.55052331 11.3130875 Oct 5 19637.92 152.61117570 6.0847258 Nov 5 19188.57 127.49813834 6.4324482 Dec 5 19679.51 142.63233452 -29.5089195 Jan 6 20692.36 178.00379600 137.6355003 Feb 6 23511.88 290.65079585 83.1214254 Mar 6 22978.64 256.05565664 -41.6376106 Apr 6 21866.00 198.75023837 -52.0010400 May 6 21876.55 190.86365569 51.4510411 Jun 6 21835.88 181.15395152 -58.8791520 Jul 6 21328.71 152.27095946 54.2924602 Aug 6 21495.02 152.86029132 -28.0157324 Sep 6 22048.91 169.70532948 3.0904465 Oct 6 22638.27 187.34170895 41.7261454 Nov 6 24266.33 247.92310242 53.6716703 Dec 6 25011.12 268.78120287 -34.1175691 Jan 7 25292.44 269.30165988 -88.4391434 > m$resid Jan Feb Mar Apr May Jun 1 0.00000000 -0.22716928 0.50969215 0.36403938 -0.81802731 0.14120092 2 -0.42699020 0.10699943 0.13432606 0.12439147 0.09004142 0.91603084 3 0.89322112 0.48365443 -0.23061570 1.40503050 1.34170138 -2.94264870 4 0.19626574 0.94977624 -0.71042410 0.15504205 -0.52218444 -0.34914394 5 2.06047619 0.88788444 -0.25667852 -2.50403529 -0.50307585 -0.14811009 6 1.28769784 3.65402350 -1.17932711 -1.95822450 -0.26927785 -0.33128372 7 0.01842522 Jul Aug Sep Oct Nov Dec 1 0.09245250 0.22595307 0.13100903 0.23282733 0.10682892 -0.34827184 2 -0.02440536 0.11474460 0.70096877 0.82345037 0.57828284 1.02080643 3 1.04172979 -0.45164618 -1.18038819 -0.41395177 0.86998621 -0.61426836 4 -0.37689846 0.04386654 0.96939747 0.72526113 0.73780276 -0.07552408 5 0.84557493 -1.67560154 0.35869413 1.51384177 -0.86162799 0.51991401 6 -0.98486049 0.02008425 0.57378497 0.60044165 2.06149034 0.71054953 7 > 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/1ez0e1292963909.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/27q0h1292963909.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/37q0h1292963909.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/40hhk1292963909.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/50hhk1292963909.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/6erft1292963909.tab") > > try(system("convert tmp/1ez0e1292963909.ps tmp/1ez0e1292963909.png",intern=TRUE)) character(0) > try(system("convert tmp/27q0h1292963909.ps tmp/27q0h1292963909.png",intern=TRUE)) character(0) > try(system("convert tmp/37q0h1292963909.ps tmp/37q0h1292963909.png",intern=TRUE)) character(0) > try(system("convert tmp/40hhk1292963909.ps tmp/40hhk1292963909.png",intern=TRUE)) character(0) > try(system("convert tmp/50hhk1292963909.ps tmp/50hhk1292963909.png",intern=TRUE)) character(0) > > > proc.time() user system elapsed 1.627 0.826 3.545