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(20366,22782,19169,13807,29743,25591,29096,26482,22405,27044,17970,18730,19684,19785,18479,10698,31956,29506,34506,27165,26736,23691,18157,17328,18205,20995,17382,9367,31124,26551,30651,25859,25100,25778,20418,18688,20424,24776,19814,12738,31566,30111,30019,31934,25826,26835,20205,17789,20520,22518,15572,11509,25447,24090,27786,26195,20516,22759,19028,16971) > 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 2529115 2627841 0 > m$fitted level slope sea Jan 1 20366.00 0.00000 0.00000 Feb 1 21192.85 958.18029 1589.15481 Mar 1 20924.51 237.62579 -1755.50756 Apr 1 16728.11 -2217.32227 -2921.10817 May 1 23289.72 2906.12969 6453.27652 Jun 1 27163.60 3456.19723 -1572.60339 Jul 1 29875.26 3038.61256 -779.25946 Aug 1 28864.11 757.88033 -2382.10542 Sep 1 24645.28 -2053.83037 -2240.27814 Oct 1 24815.49 -796.70608 2228.50550 Nov 1 20595.29 -2731.40230 -2625.29292 Dec 1 18041.78 -2630.90635 688.22256 Jan 2 17940.77 -1209.66564 1743.23420 Feb 2 16975.88 -1071.41298 2809.11943 Mar 2 17684.13 -79.20622 794.87019 Apr 2 17854.52 56.30387 -7156.51620 May 2 23169.61 2934.10363 8786.38900 Jun 2 29556.18 4836.07240 -50.17997 Jul 2 34022.21 4633.85362 483.79440 Aug 2 31726.92 863.18658 -4561.92322 Sep 2 29995.59 -551.22699 -3259.59184 Oct 2 23177.64 -3973.12434 513.36035 Nov 2 19625.17 -3743.46335 -1468.17331 Dec 2 16921.75 -3176.10879 406.24833 Jan 3 15389.44 -2278.77828 2815.56239 Feb 3 16663.23 -336.70669 4331.76668 Mar 3 17277.27 180.42962 104.73272 Apr 3 19048.07 1038.98436 -9681.07432 May 3 23112.43 2674.04732 8011.57215 Jun 3 26436.67 3027.02459 114.32964 Jul 3 27994.97 2230.62477 2656.02603 Aug 3 28999.38 1568.08966 -3140.38234 Sep 3 27138.73 -282.32475 -2038.73141 Oct 3 25488.77 -1020.78188 289.22756 Nov 3 22910.09 -1862.35775 -2492.08624 Dec 3 19653.32 -2616.00550 -965.31810 Jan 4 18146.31 -2016.04642 2277.69021 Feb 4 19263.46 -321.40537 5512.53761 Mar 4 20470.38 502.62619 -656.38201 Apr 4 23113.54 1652.88333 -10375.53524 May 4 24521.75 1521.34824 7044.25106 Jun 4 28219.11 2694.19582 1891.88701 Jul 4 28457.00 1369.86227 1561.99604 Aug 4 31873.74 2471.01458 60.25541 Sep 4 29872.06 69.71035 -4046.05609 Oct 4 26574.90 -1736.95553 260.09812 Nov 4 22654.28 -2909.39881 -2449.27738 Dec 4 19331.41 -3131.62734 -1542.41000 Jan 5 18392.80 -1952.09548 2127.20396 Feb 5 17551.70 -1354.97375 4966.29836 Mar 5 16924.79 -964.57460 -1352.79035 Apr 5 19977.37 1185.26905 -8468.37446 May 5 20501.25 831.20032 4945.74774 Jun 5 21503.69 923.02599 2586.31415 Jul 5 25315.73 2473.13875 2470.27369 Aug 5 25492.25 1242.65899 702.75349 Sep 5 23839.26 -305.71389 -3323.26081 Oct 5 21787.94 -1238.29267 971.06387 Nov 5 20894.18 -1054.13716 -1866.18209 Dec 5 19616.59 -1173.70450 -2645.59356 > m$resid Jan Feb Mar Apr May Jun 1 0.00000000 0.98896720 -0.43337536 -1.72477544 3.30016681 0.34092376 2 0.89809425 0.08842433 0.61885561 0.08543735 1.84118318 1.19792109 3 0.56595987 1.22336395 0.32416424 0.53979586 1.03437660 0.22272200 4 0.37778735 1.06524752 0.51719185 0.72330051 -0.08296574 0.73948144 5 0.74223943 0.37518910 0.24515726 1.35197119 -0.22309598 0.05785954 Jul Aug Sep Oct Nov Dec 1 -0.26247166 -1.43811313 -1.76855918 0.79036271 -1.21659674 0.06319442 2 -0.12652541 -2.36958845 -0.89016952 -2.15211395 0.14443276 0.35710662 3 -0.49994811 -0.41577601 -1.16318784 -0.46458320 -0.52967021 -0.47458871 4 -0.83258960 0.69128107 -1.50860829 -1.13655259 -0.73825236 -0.13995657 5 0.97495328 -0.77280228 -0.97267419 -0.58659192 0.11596051 -0.07529578 > 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/1is0r1261129592.ps",horizontal=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/2o1ed1261129592.ps",horizontal=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/3n8uh1261129592.ps",horizontal=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/48t8d1261129592.ps",horizontal=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/5aj011261129592.ps",horizontal=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/668pd1261129592.tab") > try(system("convert tmp/1is0r1261129592.ps tmp/1is0r1261129592.png",intern=TRUE)) character(0) > try(system("convert tmp/2o1ed1261129592.ps tmp/2o1ed1261129592.png",intern=TRUE)) character(0) > try(system("convert tmp/3n8uh1261129592.ps tmp/3n8uh1261129592.png",intern=TRUE)) character(0) > try(system("convert tmp/48t8d1261129592.ps tmp/48t8d1261129592.png",intern=TRUE)) character(0) > try(system("convert tmp/5aj011261129592.ps tmp/5aj011261129592.png",intern=TRUE)) character(0) > > > proc.time() user system elapsed 1.252 0.813 1.503