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(2849.27,2921.44,2981.85,3080.58,3106.22,3119.31,3061.26,3097.31,3161.69,3257.16,3277.01,3295.32,3363.99,3494.17,3667.03,3813.06,3917.96,3895.51,3801.06,3570.12,3701.61,3862.27,3970.1,4138.52,4199.75,4290.89,4443.91,4502.64,4356.98,4591.27,4696.96,4621.4,4562.84,4202.52,4296.49,4435.23,4105.18,4116.68,3844.49,3720.98,3674.4,3857.62,3801.06,3504.37,3032.6,3047.03,2962.34,2197.82,2014.45,1862.83,1905.41,1810.99,1670.07,1864.44,2052.02,2029.6,2070.83,2293.41,2443.27,2513.17) > 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 23691.570 1510.903 0.000 0.000 > m$fitted level slope sea Jan 1 2849.270 0.000000 0.000000 Feb 1 2917.763 5.137906 3.676866 Mar 1 2978.128 11.232362 3.722102 Apr 1 3076.797 24.202822 3.783097 May 1 3102.436 24.454221 3.783923 Jun 1 3115.531 22.264318 3.778647 Jul 1 3057.511 5.870326 3.748967 Aug 1 3093.552 12.247057 3.757763 Sep 1 3157.920 23.493537 3.769680 Oct 1 3253.377 39.216448 3.782539 Nov 1 3273.230 34.953854 3.779841 Dec 1 3291.542 31.273756 3.778034 Jan 2 3381.427 43.738211 -17.436953 Feb 2 3491.102 57.851864 3.067995 Mar 2 3664.062 83.464393 2.968218 Apr 2 3810.134 97.396664 2.926026 May 2 3915.038 99.067439 2.922092 Jun 2 3892.538 72.009896 2.971609 Jul 2 3798.036 34.946001 3.024336 Aug 2 3567.030 -24.255387 3.089805 Sep 2 3698.550 10.422176 3.059995 Oct 2 3859.232 43.873310 3.037641 Nov 2 3967.070 58.113521 3.030244 Dec 2 4135.500 82.673658 3.020327 Jan 3 4230.564 85.391652 -30.814057 Feb 3 4288.485 79.434206 2.404596 Mar 3 4441.550 95.839846 2.360016 Apr 3 4500.263 87.570663 2.377489 May 3 4354.517 35.616218 2.462843 Jun 3 4588.864 79.864375 2.406329 Jul 3 4694.559 85.615538 2.400619 Aug 3 4618.972 49.725506 2.428321 Sep 3 4560.397 25.613860 2.442788 Oct 3 4200.037 -60.318832 2.482869 Nov 3 4294.020 -25.965143 2.470413 Dec 3 4432.770 10.707486 2.460077 Jan 4 4167.371 -50.181903 -62.191487 Feb 4 4111.091 -51.513580 5.588812 Mar 4 3838.802 -100.696153 5.688161 Apr 4 3715.284 -105.778743 5.696143 May 4 3668.720 -92.593231 5.680043 Jun 4 3851.998 -31.169249 5.621739 Jul 4 3795.434 -36.823319 5.625911 Aug 4 3498.711 -94.688020 5.659104 Sep 4 3026.903 -178.650508 5.696543 Oct 4 3041.348 -135.659292 5.681641 Nov 4 2956.661 -124.310623 5.678583 Dec 4 2192.112 -266.856731 5.708441 Jan 5 2064.115 -236.158080 -49.664980 Feb 5 1858.061 -229.557601 4.769017 Mar 5 1900.738 -168.918027 4.671728 Apr 5 1806.339 -152.322988 4.651028 May 5 1665.421 -149.783354 4.648565 Jun 5 1859.849 -73.143349 4.590787 Jul 5 2047.463 -15.086448 4.556762 Aug 5 2025.042 -16.719388 4.557506 Sep 5 2066.277 -3.816287 4.552936 Oct 5 2288.871 46.592705 4.539059 Nov 5 2438.736 69.585808 4.534138 Dec 5 2508.636 69.655764 4.534127 > m$resid Jan Feb Mar Apr May Jun 1 0.000000000 0.264174668 0.338719245 0.524636601 0.008481258 -0.066337608 2 0.391778289 0.317896052 0.658757312 0.358375196 0.042979359 -0.696059102 3 0.076849005 -0.142027442 0.421427347 -0.212544699 -1.335876533 1.137975849 4 -1.668936108 -0.032468156 -1.263878498 -0.130668861 0.339078839 1.579836454 5 0.828389947 0.162901086 1.558656848 0.426703256 0.065314781 1.971295881 Jul Aug Sep Oct Nov Dec 1 -0.465499328 0.174130836 0.299950075 0.413402697 -0.111114812 -0.095432312 2 -0.953494777 -1.523017490 0.892123635 0.860576335 0.366350051 0.631846223 3 0.147928026 -0.923215445 -0.620264419 -2.210655129 0.883778750 0.943445458 4 -0.145438009 -1.488525210 -2.159948123 1.105979691 0.291956419 -3.667177163 5 1.493426910 -0.042006936 0.331938029 1.296816087 0.591524113 0.001799721 > 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/16w131259949408.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/21eyj1259949408.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/36xlo1259949408.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/4hdwd1259949408.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/5h5xs1259949408.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/6ayj81259949408.tab") > system("convert tmp/16w131259949408.ps tmp/16w131259949408.png") > system("convert tmp/21eyj1259949408.ps tmp/21eyj1259949408.png") > system("convert tmp/36xlo1259949408.ps tmp/36xlo1259949408.png") > system("convert tmp/4hdwd1259949408.ps tmp/4hdwd1259949408.png") > system("convert tmp/5h5xs1259949408.ps tmp/5h5xs1259949408.png") > > > proc.time() user system elapsed 1.431 0.819 1.805