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(19915,19843,19761,20858,21968,23061,22661,22269,21857,21568,21274,20987,19683,19381,19071,20772,22485,24181,23479,22782,22067,21489,20903,20330,19736,19483,19242,20334,21423,22523,21986,21462,20908,20575,20237,19904,19610,19251,18941,20450,21946,23409,22741,22069,21539,21189,20960,20704,19697,19598,19456,20316,21083,22158,21469,20892,20578,20233,19947,20049) > 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.000 362973.077 3396.667 0.000 > m$fitted level slope sea Jan 1 19915.00 0.00000 0.0000000 Feb 1 19843.99 -71.08786 -0.9943179 Mar 1 19760.73 -82.79603 0.2690167 Apr 1 20822.08 1008.51469 35.9214797 May 1 21973.67 1145.05558 -5.6653171 Jun 1 23064.75 1093.53921 -3.7500111 Jul 1 22707.95 -290.80900 -46.9502450 Aug 1 22262.82 -438.11154 6.1842030 Sep 1 21855.58 -408.64357 1.4228865 Oct 1 21565.30 -295.66471 2.7013828 Nov 1 21275.63 -289.94319 -1.6281968 Dec 1 20987.82 -287.90180 -0.8237164 Jan 2 19713.57 -1229.09329 -30.5671244 Feb 2 19310.59 -440.18827 70.4115254 Mar 2 19115.79 -211.95664 -44.7896344 Apr 2 20656.38 1417.39023 115.6227341 May 2 22522.44 1834.19200 -37.4357335 Jun 2 24141.51 1634.28478 39.4906967 Jul 2 23624.43 -364.96891 -145.4298295 Aug 2 22775.85 -814.37350 6.1471578 Sep 2 22056.27 -726.28766 10.7305202 Oct 2 21478.67 -588.11838 10.3331100 Nov 2 20933.57 -548.13915 -30.5700213 Dec 2 20275.40 -650.38813 54.6024161 Jan 3 19789.53 -497.55170 -53.5287823 Feb 3 19358.70 -435.51783 124.3048773 Mar 3 19360.82 -34.35954 -118.8158324 Apr 3 20212.04 779.02827 121.9614160 May 3 21503.32 1248.84228 -80.3169350 Jun 3 22377.38 904.97855 145.6249821 Jul 3 22155.86 -128.62536 -169.8645339 Aug 3 21477.50 -633.02258 -15.5015573 Sep 3 20895.71 -586.01573 12.2880890 Oct 3 20556.16 -359.88405 18.8440521 Nov 3 20254.93 -306.06386 -17.9294447 Dec 3 19868.31 -379.96974 35.6895199 Jan 4 19643.49 -237.65630 -33.4889262 Feb 4 19137.84 -483.58349 113.1555129 Mar 4 19081.91 -94.15839 -140.9069837 Apr 4 20312.20 1114.45402 137.8031822 May 4 22015.43 1651.00213 -69.4320284 Jun 4 23226.23 1249.74236 182.7729538 Jul 4 22952.31 -139.24447 -211.3107253 Aug 4 22111.02 -779.23172 -42.0188323 Sep 4 21523.76 -604.22883 15.2367403 Oct 4 21158.67 -386.23592 30.3321942 Nov 4 20961.77 -213.63958 -1.7651677 Dec 4 20701.53 -256.11333 2.4727916 Jan 5 19728.63 -909.41822 -31.6285921 Feb 5 19424.11 -358.09496 173.8913392 Mar 5 19651.48 173.13619 -195.4821300 Apr 5 20235.67 546.74631 80.3260710 May 5 21160.34 889.87803 -77.3438615 Jun 5 21880.09 735.39115 277.9115110 Jul 5 21665.69 -127.16339 -196.6914021 Aug 5 20968.74 -644.62138 -76.7377564 Sep 5 20556.70 -433.40393 21.3020191 Oct 5 20221.43 -344.27839 11.5662702 Nov 5 19964.38 -265.06936 -17.3764038 Dec 5 19945.15 -41.84546 103.8531694 > m$resid Jan Feb Mar Apr May Jun 1 0.000000000 -0.118782435 -0.019797914 1.809303118 0.226651829 -0.085508161 2 -1.562887975 1.315442851 0.382814742 2.701008170 0.691915353 -0.331810702 3 0.253817092 0.103160025 0.669222744 1.349313679 0.779813141 -0.570766242 4 0.236367687 -0.408410334 0.648006002 2.006024150 0.890443160 -0.666052638 5 -1.085063798 0.915106240 0.882959553 0.620273995 0.569405804 -0.256435029 Jul Aug Sep Oct Nov Dec 1 -2.297778744 -0.244496715 0.048911736 0.187525347 0.009496724 0.003388352 2 -3.318418682 -0.745933808 0.146207240 0.229337088 0.066358632 -0.169715975 3 -1.715601329 -0.837212479 0.078023237 0.375339452 0.089332474 -0.122671610 4 -2.305467152 -1.062268896 0.290474498 0.361830647 0.286481277 -0.070500453 5 -1.431686033 -0.858890526 0.350584444 0.147933217 0.131473536 0.370530374 > 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/1rpu21260973133.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/28d6o1260973133.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/3r8y91260973133.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/4n6b41260973133.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/5ghj71260973133.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/62rqq1260973133.tab") > > try(system("convert tmp/1rpu21260973133.ps tmp/1rpu21260973133.png",intern=TRUE)) character(0) > try(system("convert tmp/28d6o1260973133.ps tmp/28d6o1260973133.png",intern=TRUE)) character(0) > try(system("convert tmp/3r8y91260973133.ps tmp/3r8y91260973133.png",intern=TRUE)) character(0) > try(system("convert tmp/4n6b41260973133.ps tmp/4n6b41260973133.png",intern=TRUE)) character(0) > try(system("convert tmp/5ghj71260973133.ps tmp/5ghj71260973133.png",intern=TRUE)) character(0) > > > proc.time() user system elapsed 1.345 0.842 3.399