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(11100,8962,9173,8738,8459,8078,8411,8291,7810,8616,8312,9692,9911,8915,9452,9112,8472,8230,8384,8625,8221,8649,8625,10443,10357,8586,8892,8329,8101,7922,8120,7838,7735,8406,8209,9451,10041,9411,10405,8467,8464,8102,7627,7513,7510,8291,8064,9383,9706,8579,9474,8318,8213,8059,9111,7708,7680,8014,8007,8718,9486,9113,9025,8476,7952,7759,7835,7600,7651,8319,8812,8630) > 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 198812.06935 80.40304 16369.53784 0.00000 > m$fitted level slope sea Jan 1 11100.000 0.00000 0.000000 Feb 1 9151.729 -104.28120 -189.728988 Mar 1 9190.941 -101.30990 -17.941392 Apr 1 8888.553 -102.67484 -150.553377 May 1 8586.609 -103.88692 -127.608970 Jun 1 8219.341 -105.73308 -141.341206 Jul 1 8434.167 -103.31493 -23.167412 Aug 1 8401.905 -102.75402 -110.904936 Sep 1 7980.192 -105.37678 -170.192131 Oct 1 8562.730 -99.49347 53.269636 Nov 1 8460.451 -99.51821 -148.450871 Dec 1 9557.148 -88.51455 134.851564 Jan 2 9113.967 -76.79942 797.033276 Feb 2 9109.202 -75.22518 -194.202029 Mar 2 9337.537 -69.13439 114.463307 Apr 2 9234.076 -69.53022 -122.075690 May 2 8660.213 -73.98935 -188.213181 Jun 2 8404.867 -75.62865 -174.866814 Jul 2 8361.546 -75.31758 22.454198 Aug 2 8564.600 -72.51941 60.400419 Sep 2 8483.619 -72.60752 -262.619347 Oct 2 8519.855 -71.42356 129.145168 Nov 2 8846.302 -67.04416 -221.302009 Dec 2 9950.458 -58.06694 492.542494 Jan 3 9755.945 -57.61577 601.054977 Feb 3 9043.113 -68.14919 -457.113402 Mar 3 8787.633 -71.79211 104.366568 Apr 3 8415.528 -76.39118 -86.527709 May 3 8250.084 -77.51468 -149.083769 Jun 3 8113.146 -78.23480 -191.145889 Jul 3 8123.307 -77.13476 -3.306612 Aug 3 7832.570 -79.88440 5.430327 Sep 3 7934.187 -77.47330 -199.186877 Oct 3 8262.374 -71.96526 143.625924 Nov 3 8611.015 -66.49924 -402.015291 Dec 3 8868.111 -63.21895 582.889081 Jan 4 9171.825 -60.59658 869.175133 Feb 4 9683.799 -51.56756 -272.798903 Mar 4 10131.993 -41.98669 273.006650 Apr 4 8930.001 -62.27646 -463.001384 May 4 8618.654 -66.13527 -154.653554 Jun 4 8339.830 -69.27110 -237.829814 Jul 4 7714.381 -77.50524 -87.381069 Aug 4 7529.725 -79.12653 -16.724987 Sep 4 7700.836 -75.25963 -190.836048 Oct 4 8099.916 -67.90982 191.083946 Nov 4 8441.922 -61.87519 -377.922003 Dec 4 8800.756 -56.43320 582.244325 Jan 5 8962.063 -53.74298 743.936738 Feb 5 8958.933 -52.89831 -379.932949 Mar 5 8936.913 -52.30553 537.086871 Apr 5 8771.651 -54.41457 -453.651225 May 5 8400.380 -59.92170 -187.379998 Jun 5 8182.956 -62.54710 -123.956045 Jul 5 8905.580 -49.53749 205.419862 Aug 5 8062.463 -62.82988 -354.462824 Sep 5 7956.007 -63.56828 -276.006804 Oct 5 7911.391 -63.24957 102.608977 Nov 5 8334.512 -55.39514 -327.511557 Dec 5 8228.535 -56.16169 489.465016 Jan 6 8617.069 -49.36920 868.930808 Feb 6 9302.619 -36.39437 -189.618862 Mar 6 8678.008 -47.77494 346.991638 Apr 6 8788.094 -44.72823 -312.094084 May 6 8285.387 -53.21370 -333.386920 Jun 6 8076.973 -56.00111 -317.972624 Jul 6 7565.132 -64.11034 269.868356 Aug 6 7810.255 -58.59176 -210.255042 Sep 6 7898.112 -55.97049 -247.112260 Oct 6 8212.679 -49.39984 106.321069 Nov 6 8922.789 -36.30429 -110.789234 Dec 6 8450.529 -43.59470 179.470934 > m$resid Jan Feb Mar Apr May Jun 1 0.000000000 -2.682748378 0.316778501 -0.451974963 -0.447458381 -0.591161030 2 -0.919360847 0.145083277 0.659337042 -0.076785124 -1.128177931 -0.405439422 3 -0.319682586 -1.405665994 -0.406621489 -0.668276694 -0.198626920 -0.132452644 4 0.833591910 1.250940958 1.088752465 -2.570978387 -0.554155911 -0.472992881 5 0.488438298 0.111257269 0.067514516 -0.249736023 -0.703551818 -0.349681683 6 0.991526313 1.619312133 -1.289658731 0.348560650 -1.015362181 -0.344168948 Jul Aug Sep Oct Nov Dec 1 0.719302188 0.159395276 -0.715355767 1.542456759 -0.006245305 2.680846570 2 0.072200114 0.621935553 -0.018899659 0.243132176 0.888633715 2.612693917 3 0.196954523 -0.475823843 0.404301991 0.903678996 0.936179493 0.720963365 4 -1.236400091 -0.238168978 0.556232192 1.054178990 0.910033594 0.936007505 5 1.742572240 -1.761082737 -0.096811542 0.042039321 1.077946497 -0.112371213 6 -1.010505237 0.685424121 0.324561289 0.820716059 1.681519546 -0.967323187 > 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/1ivcp1293542031.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/2ivcp1293542031.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/3mwsd1293542031.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/4mwsd1293542031.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/5mwsd1293542031.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/6ioqm1293542031.tab") > > try(system("convert tmp/1ivcp1293542031.ps tmp/1ivcp1293542031.png",intern=TRUE)) character(0) > try(system("convert tmp/2ivcp1293542031.ps tmp/2ivcp1293542031.png",intern=TRUE)) character(0) > try(system("convert tmp/3mwsd1293542031.ps tmp/3mwsd1293542031.png",intern=TRUE)) character(0) > try(system("convert tmp/4mwsd1293542031.ps tmp/4mwsd1293542031.png",intern=TRUE)) character(0) > try(system("convert tmp/5mwsd1293542031.ps tmp/5mwsd1293542031.png",intern=TRUE)) character(0) > > > proc.time() user system elapsed 1.543 0.845 3.433