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(785.8,819.3,849.4,880.4,900.1,937.2,948.9,952.6,947.3,974.2,1000.8,1032.8,1050.7,1057.3,1075.4,1118.4,1179.8,1227,1257.8,1251.5,1236.3,1170.6,1213.1,1265.5,1300.8,1348.4,1371.9,1403.3,1451.8,1474.2,1438.2,1513.6,1562.2,1546.2,1527.5,1418.7,1448.5,1492.1,1395.4,1403.7,1316.6,1274.5,1264.4,1323.9,1332.1,1250.2,1096.7,1080.8,1039.2,792,746.6,688.8,715.8,672.9,629.5,681.2,755.4,760.6,765.9,836.8,904.9) > 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 2236.9490 209.4151 0.0000 0.0000 > m$fitted level slope sea Jan 1 785.8000 0.000000 0.0000000 Feb 1 817.6093 2.676917 1.6907176 Mar 1 847.7204 6.845772 1.6795776 Apr 1 878.7283 11.608924 1.6717003 May 1 898.4303 13.431803 1.6696561 Jun 1 935.5349 19.154747 1.6651223 Jul 1 947.2338 17.282344 1.6661916 Aug 1 950.9324 13.800006 1.6676404 Sep 1 945.6309 8.848117 1.6691496 Oct 1 972.5319 13.556676 1.6680953 Nov 1 999.1325 16.970030 1.6675329 Dec 1 1031.1329 20.910308 1.6670547 Jan 2 1064.1953 23.986108 -13.4952809 Feb 2 1056.7748 16.140104 0.5251740 Mar 2 1074.8772 16.656275 0.5228219 Apr 2 1117.9005 23.586545 0.4995243 May 2 1179.3251 33.528100 0.4748753 Jun 2 1226.5317 37.121394 0.4683054 Jul 2 1257.3295 35.460284 0.4705449 Aug 2 1251.0185 24.487762 0.4814518 Sep 2 1235.8109 14.060352 0.4890940 Oct 2 1170.0996 -6.894970 0.5004174 Nov 2 1212.6048 6.082295 0.4952473 Dec 2 1265.0083 18.250983 0.4916729 Jan 3 1302.5884 23.252248 -1.7884357 Feb 3 1347.9204 28.866883 0.4796461 Mar 3 1371.4161 27.454226 0.4839121 Apr 3 1402.8184 28.491926 0.4816005 May 3 1451.3270 33.751379 0.4729603 Jun 3 1473.7234 30.768243 0.4765740 Jul 3 1437.7078 13.224024 0.4922445 Aug 3 1513.1185 29.560333 0.4814859 Sep 3 1561.7209 34.562669 0.4790570 Oct 3 1545.7162 21.278518 0.4838127 Nov 3 1527.0134 10.775230 0.4865850 Dec 3 1418.2073 -20.639759 0.4926985 Jan 4 1447.7761 -7.579106 0.7238591 Feb 4 1491.6214 5.611236 0.4786284 Mar 4 1394.8611 -21.305670 0.5388546 Apr 4 1403.1740 -13.521747 0.5260091 May 4 1316.0505 -32.860207 0.5495430 Jun 4 1273.9483 -35.288228 0.5517217 Jul 4 1263.8527 -28.669979 0.5473430 Aug 4 1323.3640 -5.504464 0.5360424 Sep 4 1331.5653 -1.903904 0.5347474 Oct 4 1249.6597 -22.920851 0.5403206 Nov 4 1096.1530 -57.226905 0.5470278 Dec 4 1080.2545 -46.369459 0.5454627 Jan 5 1044.2576 -43.664347 -5.0576322 Feb 5 793.2784 -97.108840 -1.2784379 Mar 5 747.9026 -83.508775 -1.3025992 Apr 5 690.1115 -76.750487 -1.3114540 May 5 717.1378 -49.484236 -1.3377973 Jun 5 674.2390 -47.754119 -1.3390298 Jul 5 630.8396 -46.610096 -1.3396307 Aug 5 682.5496 -20.780752 -1.3496336 Sep 5 756.7568 4.173260 -1.3567587 Oct 5 761.9568 4.443009 -1.3568154 Nov 5 767.2569 4.668159 -1.3568504 Dec 5 838.1588 22.068635 -1.3588416 Jan 6 887.9496 29.309346 16.9504229 > m$resid Jan Feb Mar Apr May Jun 1 0.00000000 0.39719223 0.53392600 0.45758259 0.15057277 0.43574192 2 0.25614210 -0.47685659 0.03556565 0.47814759 0.68640183 0.24819135 3 0.37896552 0.36046037 -0.09740109 0.07162107 0.36320375 -0.20606936 4 0.96007478 0.86551753 -1.85696181 0.53740806 -1.33569011 -0.16773865 5 0.19583242 -3.54866067 0.93857040 0.46668392 1.88344525 0.11953081 6 0.51936801 Jul Aug Sep Oct Nov Dec 1 -0.13639445 -0.24763963 -0.34756660 0.32814458 0.23696193 0.27296739 2 -0.11475842 -0.75812837 -0.72051005 -1.44801391 0.89674634 0.84088109 3 -1.21211966 1.12876620 0.34565585 -0.91794418 -0.72579489 -2.17084572 4 0.45727396 1.60067780 0.24879840 -1.45229757 -2.37061674 0.75027551 5 0.07904622 1.78477046 1.72433530 0.01864008 0.01555837 1.20241618 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/10y1f1259693618.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/24fsv1259693618.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/30pvs1259693618.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/4d7ya1259693618.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/52uys1259693618.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/6s8ym1259693618.tab") > > system("convert tmp/10y1f1259693618.ps tmp/10y1f1259693618.png") > system("convert tmp/24fsv1259693618.ps tmp/24fsv1259693618.png") > system("convert tmp/30pvs1259693618.ps tmp/30pvs1259693618.png") > system("convert tmp/4d7ya1259693618.ps tmp/4d7ya1259693618.png") > system("convert tmp/52uys1259693618.ps tmp/52uys1259693618.png") > > > proc.time() user system elapsed 1.288 0.815 1.784