R version 2.12.1 (2010-12-16) Copyright (C) 2010 The R Foundation for Statistical Computing ISBN 3-900051-07-0 Platform: i486-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(2564,2820,3508,3088,3299,2939,3320,3418,3604,3495,4163,4882,2211,3260,2992,2425,2707,3244,3965,3315,3333,3583,4021,4904,2252,2952,3573,3048,3059,2731,3563,3092,3478,3478,4308,5029,2075,3264,3308,3688,3136,2824,3644,4694,2914,3686,4358,5587,2265,3685,3754,3708,3210,3517,3905,3670,4221,4404,5086,5725) > 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 50437.58 0.00 61303.99 48227.42 > m$fitted level slope sea Jan 1 2564.000 0.000000 0.0000000 Feb 1 2688.987 5.834267 88.9692422 Mar 1 3112.632 34.712980 290.6774690 Apr 1 3164.116 35.634465 -81.3838736 May 1 3228.592 36.737271 60.1373753 Jun 1 3123.472 32.814695 -131.6657611 Jul 1 3182.071 33.394397 128.1920018 Aug 1 3288.838 34.885352 101.3342030 Sep 1 3437.235 37.095171 123.6644917 Oct 1 3489.843 37.390937 -0.7339086 Nov 1 3772.258 41.990837 297.6746345 Dec 1 4276.959 50.542703 429.2588198 Jan 2 3838.527 70.972931 -1430.6071707 Feb 2 3633.455 70.059818 -274.0145631 Mar 2 3226.037 57.255877 -83.7542774 Apr 2 2869.160 44.479688 -312.8270170 May 2 2697.463 38.674621 81.5796183 Jun 2 2935.252 42.936280 240.0544522 Jul 2 3346.337 49.286215 489.4816029 Aug 2 3416.572 49.596230 -108.9766152 Sep 2 3399.463 48.696048 -42.8172580 Oct 2 3561.910 50.107500 -19.2903737 Nov 2 3739.288 51.408151 236.4528205 Dec 2 3948.604 52.186159 898.9557874 Jan 3 3797.005 53.467886 -1470.5594115 Feb 3 3433.490 51.964894 -334.0440333 Mar 3 3367.361 50.194367 245.1262828 Apr 3 3340.614 48.644426 -267.2973141 May 3 3278.822 46.408528 -182.9075271 Jun 3 3064.744 41.722880 -244.8328869 Jul 3 3052.388 40.894524 529.3457932 Aug 3 3124.510 41.302536 -43.4190275 Sep 3 3344.027 43.304346 71.4425654 Oct 3 3498.360 44.344450 -59.4088016 Nov 3 3773.788 45.943060 452.7842617 Dec 3 3902.872 46.253754 1096.7376523 Jan 4 3717.285 45.961338 -1559.9898286 Feb 4 3625.053 45.269309 -312.7324537 Mar 4 3378.939 42.144871 28.3398841 Apr 4 3523.813 43.639639 129.6449418 May 4 3427.102 41.434178 -243.7222853 Jun 4 3288.924 38.741807 -403.5857626 Jul 4 3243.200 37.611434 429.9661745 Aug 4 3860.342 44.302813 632.0051393 Sep 4 3651.710 41.847971 -649.2946670 Oct 4 3722.972 42.077978 -47.2797692 Nov 4 3818.140 42.388919 521.2229509 Dec 4 4037.489 43.120169 1487.3095611 Jan 5 3964.707 42.700565 -1659.0029263 Feb 5 3926.028 42.244370 -212.7038286 Mar 5 3867.258 41.356338 -78.5803482 Apr 5 3723.464 39.246162 47.5390729 May 5 3584.911 37.003215 -314.3901854 Jun 5 3725.794 38.304136 -244.3791810 Jul 5 3773.649 38.414581 128.0547021 Aug 5 3481.483 35.054063 303.3452710 Sep 5 3980.170 39.038226 79.1184357 Oct 5 4277.739 40.847115 35.8668073 Nov 5 4478.205 41.731606 551.9111262 Dec 5 4408.675 41.225260 1355.2917291 > m$resid Jan Feb Mar Apr May Jun 1 0.00000000 0.47241613 1.36695518 0.06528600 0.12236502 -0.61906544 2 -2.45003234 -1.21720302 -1.90869055 -1.68546500 -0.91629489 0.86548830 3 -0.92746170 -1.83962288 -0.49917047 -0.32316377 -0.47137691 -1.13099169 4 -1.03289927 -0.60749292 -1.25563289 0.43949089 -0.60410228 -0.78122121 5 -0.51306475 -0.35746891 -0.43907505 -0.80061607 -0.77058581 0.45305192 Jul Aug Sep Oct Nov Dec 1 0.11347024 0.32362990 0.50091345 0.06845434 1.08128410 2.04204759 2 1.61850667 0.09251446 -0.29506547 0.50339361 0.56321503 0.70063612 3 -0.23738775 0.13785605 0.78868721 0.49170923 1.02364503 0.36901535 4 -0.37060770 2.55656903 -1.11878545 0.13022641 0.23514631 0.78419533 5 0.04191689 -1.45753983 2.04930944 1.14397243 0.70654283 -0.49237930 > 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/rcomp/tmp/1ixuw1322499080.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/rcomp/tmp/2ysmd1322499080.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/rcomp/tmp/3p8751322499080.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/rcomp/tmp/41mnx1322499080.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/rcomp/tmp/5x1h61322499080.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/rcomp/createtable file can be downloaded at http://www.wessa.net/cretab > load(file="/var/www/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/rcomp/tmp/6ujxc1322499080.tab") > > try(system("convert tmp/1ixuw1322499080.ps tmp/1ixuw1322499080.png",intern=TRUE)) character(0) > try(system("convert tmp/2ysmd1322499080.ps tmp/2ysmd1322499080.png",intern=TRUE)) character(0) > try(system("convert tmp/3p8751322499080.ps tmp/3p8751322499080.png",intern=TRUE)) character(0) > try(system("convert tmp/41mnx1322499080.ps tmp/41mnx1322499080.png",intern=TRUE)) character(0) > try(system("convert tmp/5x1h61322499080.ps tmp/5x1h61322499080.png",intern=TRUE)) character(0) > > > proc.time() user system elapsed 2.232 0.388 2.627