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(462,455,461,461,463,462,456,455,456,472,472,471,465,459,465,468,467,463,460,462,461,476,476,471,453,443,442,444,438,427,424,416,406,431,434,418,412,404,409,412,406,398,397,385,390,413,413,401,397,397,409,419,424,428,430,424,433,456,459,446,441) > 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 4.067326 7.484846 0.000000 > m$fitted level slope sea Jan 1 462.0000 0.0000000 0.0000000 Feb 1 460.6461 -1.8205545 -5.6460622 Mar 1 457.6985 -2.4639820 3.3014735 Apr 1 458.0381 -1.0821704 2.9618507 May 1 460.5189 0.8048590 2.4811295 Jun 1 462.1862 1.2579157 -0.1862395 Jul 1 459.1136 -0.9566004 -3.1135667 Aug 1 455.4810 -2.3214055 -0.4810121 Sep 1 454.1651 -1.8057974 1.8349459 Oct 1 463.5395 3.9416403 8.4604796 Nov 1 471.6371 6.0785457 0.3629261 Dec 1 474.3434 4.3451159 -3.3434248 Jan 2 470.7662 0.3010715 -5.7661535 Feb 2 466.5588 -2.0136841 -7.5587957 Mar 2 463.6952 -2.4472527 1.3047957 Apr 2 464.4973 -0.8342239 3.5027184 May 2 464.1018 -0.6166543 2.8982266 Jun 2 461.2253 -1.7510927 1.7746630 Jul 2 460.1170 -1.4293685 -0.1169725 Aug 2 462.2434 0.3362437 -0.2433837 Sep 2 464.9973 1.5346347 -3.9972580 Oct 2 468.6790 2.6007659 7.3210092 Nov 2 472.6034 3.2585231 3.3966234 Dec 2 472.5560 1.6158532 -1.5560424 Jan 3 463.6286 -3.6287485 -10.6285926 Feb 3 453.3213 -6.9554031 -10.3213315 Mar 3 443.1222 -8.5648874 -1.1221705 Apr 3 437.1842 -7.2711362 6.8157748 May 3 432.2548 -6.1195139 5.7451533 Jun 3 426.0952 -6.1393230 0.9047607 Jul 3 423.5868 -4.3421496 0.4132136 Aug 3 419.3424 -4.2939126 -3.3424184 Sep 3 413.5629 -5.0240769 -7.5628562 Oct 3 417.7209 -0.5167794 13.2791474 Nov 3 424.4512 3.0430190 9.5487599 Dec 3 420.0581 -0.6148303 -2.0581234 Jan 4 417.8986 -1.3756921 -5.8986445 Feb 4 414.3278 -2.4566451 -10.3278367 Mar 4 411.4626 -2.6572523 -2.4625726 Apr 4 407.1876 -3.4485403 4.8123673 May 4 401.5041 -4.5410977 4.4958940 Jun 4 397.3780 -4.3377243 0.6219815 Jul 4 394.9291 -3.4105760 2.0709393 Aug 4 390.5923 -3.8645527 -5.5922735 Sep 4 396.4728 0.8970245 -6.4728043 Oct 4 401.6171 2.9681317 11.3829119 Nov 4 402.1735 1.7918508 10.8265106 Dec 4 402.5420 1.0965514 -1.5420264 Jan 5 402.4384 0.5095376 -5.4384084 Feb 5 405.0231 1.5239011 -8.0231139 Mar 5 408.7981 2.6215876 0.2018935 Apr 5 412.3442 3.0713206 6.6558179 May 5 417.7397 4.2015992 6.2603468 Jun 5 425.7042 6.0345452 2.2958313 Jul 5 430.7111 5.5333585 -0.7111370 Aug 5 435.5452 5.1925750 -11.5452246 Sep 5 440.8979 5.2704068 -7.8978785 Oct 5 444.2838 4.3558411 11.7161705 Nov 5 447.4227 3.7652512 11.5772752 Dec 5 448.6159 2.5153833 -2.6158929 Jan 6 448.7623 1.3630290 -7.7622881 > m$resid Jan Feb Mar Apr May Jun 1 0.000000000 -1.880096994 -0.300541729 0.760951406 1.000452793 0.221881565 2 -2.016580890 -1.167157659 -0.213843572 0.795602610 0.109751572 -0.568265462 3 -2.608721575 -1.651977482 -0.795516256 0.640253123 0.573898753 -0.009883718 4 -0.377953359 -0.535701300 -0.099239104 -0.392024163 -0.543213720 0.101239213 5 -0.291453933 0.502570518 0.543311935 0.222886485 0.561483646 0.911356616 6 -0.571992137 Jul Aug Sep Oct Nov Dec 1 -1.087723385 -0.678839082 0.256269280 2.849927725 1.059398679 -0.859560870 2 0.158858310 0.871853395 0.594217921 0.529129020 0.326427003 -0.815926388 3 0.891808474 0.023849387 -0.361329735 2.235704514 1.768650481 -1.818499944 4 0.460372966 -0.224747196 2.356281532 1.026824143 -0.584405083 -0.345655369 5 -0.248819626 -0.168814524 0.038530765 -0.453358405 -0.293295965 -0.621076433 6 > 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/16htf1259943374.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/2qzsm1259943374.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/3q29y1259943374.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/4ude91259943374.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/568lb1259943374.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/6tmaj1259943374.tab") > system("convert tmp/16htf1259943374.ps tmp/16htf1259943374.png") > system("convert tmp/2qzsm1259943374.ps tmp/2qzsm1259943374.png") > system("convert tmp/3q29y1259943374.ps tmp/3q29y1259943374.png") > system("convert tmp/4ude91259943374.ps tmp/4ude91259943374.png") > system("convert tmp/568lb1259943374.ps tmp/568lb1259943374.png") > > > proc.time() user system elapsed 1.170 0.817 1.504