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(1203.6,1180.59,1156.85,1191.5,1191.33,1234.18,1220.33,1228.81,1207.01,1249.48,1248.29,1280.08,1280.66,1294.87,1310.61,1270.09,1270.2,1276.66,1303.82,1335.85,1377.94,1400.63,1418.3,1438.24,1406.82,1420.86,1482.37,1530.62,1503.35,1455.27,1473.99,1526.75,1549.38,1481.14,1468.36,1378.55,1330.63,1322.7,1385.59,1400.38,1280,1267.38,1282.83,1166.36,968.75,896.24,903.25,825.88,735.09,797.87,872.81,919.14,919.32,987.48,1020.62,1057.08,1036.19,1095.63,1115.1,1073.87,1104.49,1169.43,1186.69,1089.41,1030.71,1101.6,1049.33,1141.2,1183.26,1180.55) > 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 2093.78089 85.41624 11.92216 0.00000 > m$fitted level slope sea Jan 1 1203.6000 0.00000000 0.00000000 Feb 1 1181.8305 -1.47498088 -1.24054197 Mar 1 1158.4028 -3.12916108 -1.55276493 Apr 1 1192.0142 0.66940129 -0.51424020 May 1 1192.7061 0.67223517 -1.37611935 Jun 1 1234.7029 6.58002257 -0.52291656 Jul 1 1221.9729 3.58378469 -1.64292549 Aug 1 1229.8119 4.28068176 -1.00194943 Sep 1 1208.5773 -0.05031462 -1.56725073 Oct 1 1249.9086 7.14287249 -0.42859444 Nov 1 1249.7402 5.85172229 -1.45023735 Dec 1 1280.8229 10.35466723 -0.74290190 Jan 2 1275.2106 7.64671774 5.44936453 Feb 2 1294.5475 9.71838082 0.32250061 Mar 2 1311.5566 11.03536908 -0.94659085 Apr 2 1270.9029 1.67084498 -0.81294480 May 2 1270.9518 1.37662636 -0.75184494 Jun 2 1275.8387 2.01395506 0.82126617 Jul 2 1304.0378 6.77107289 -0.21784468 Aug 2 1334.8945 11.14856801 0.95554396 Sep 2 1378.5565 17.05950647 -0.61650847 Oct 2 1400.0701 17.86942125 0.55988538 Nov 2 1419.2902 18.11506418 -0.99017890 Dec 2 1437.8554 18.19690306 0.38461695 Jan 3 1413.6208 10.58329167 -6.80079459 Feb 3 1420.5955 9.93827699 0.26450696 Mar 3 1479.0505 18.73755282 3.31951041 Apr 3 1529.9120 24.57501853 0.70801997 May 3 1505.4338 15.66348082 -2.08376495 Jun 3 1456.1031 3.85712254 -0.83312924 Jul 3 1472.5825 6.14989163 1.40750328 Aug 3 1524.2785 14.42295407 2.47146555 Sep 3 1549.4384 16.37319932 -0.05836503 Oct 3 1483.1274 1.35429076 -1.98738189 Nov 3 1468.6042 -1.52994287 -0.24415017 Dec 3 1379.4097 -17.44342796 -0.85970929 Jan 4 1340.0929 -21.38806185 -9.46292394 Feb 4 1323.5754 -20.51383412 -0.87543920 Mar 4 1380.0651 -6.56943057 5.52485353 Apr 4 1397.7612 -2.16523547 2.61882758 May 4 1284.2396 -22.37262653 -4.23960953 Jun 4 1266.0799 -21.60820437 1.30007884 Jul 4 1280.9113 -14.99663868 1.91865127 Aug 4 1169.5417 -32.48208182 -3.18167045 Sep 4 971.7300 -62.47840002 -2.98002561 Oct 4 896.4357 -64.80368814 -0.19567776 Nov 4 896.3268 -53.06501518 6.92320581 Dec 4 827.1176 -55.99164770 -1.23763586 Jan 5 748.6173 -60.05730855 -13.52732025 Feb 5 797.7344 -40.41292405 0.13556812 Mar 5 864.0051 -21.11512384 8.80492509 Apr 5 908.0161 -9.30619823 11.12389358 May 5 922.2651 -5.03549660 -2.94505500 Jun 5 985.3269 7.30965135 2.15305850 Jul 5 1014.8191 11.33088894 5.80092345 Aug 5 1050.5177 15.74819214 6.56226703 Sep 5 1038.4752 10.71040554 -2.28515658 Oct 5 1095.2195 19.05548645 0.41046163 Nov 5 1106.9115 17.72060588 8.18848023 Dec 5 1075.5213 8.82681568 -1.65133933 Jan 6 1125.3589 16.23948105 -20.86890701 Feb 6 1172.9619 21.89113029 -3.53189607 Mar 6 1181.9975 19.56747694 4.69254037 Apr 6 1084.9747 -1.55693245 4.43534910 May 6 1037.8976 -9.80451387 -7.18762563 Jun 6 1095.7863 2.45916866 5.81372737 Jul 6 1048.3426 -6.58118111 0.98736514 Aug 6 1127.5586 8.96150787 13.64137792 Sep 6 1186.4261 18.00227419 -3.16609458 Oct 6 1182.2248 13.97981613 -1.67475562 > m$resid Jan Feb Mar Apr May 1 0.0000000000 -0.2840241986 -0.4623239297 0.7606168748 0.0004592858 2 -0.3683275145 0.1962934900 0.1437514763 -1.0170269011 -0.0319219052 3 -0.9014420636 -0.0649279144 0.9544650661 0.6308437656 -0.9636626370 4 -0.4511187952 0.0902442464 1.5125657843 0.4760599883 -2.1851124994 5 -0.4568623179 2.0560457506 2.0925690922 1.2768553530 0.4618054374 6 0.8241074893 0.5967151043 -0.2518463048 -2.2848325346 -0.8918332794 Jun Jul Aug Sep Oct 1 0.8358296127 -0.3876308784 0.0849579865 -0.5074818252 0.8208605957 2 0.0690879404 0.5153637878 0.4740413563 0.6399222299 0.0876659465 3 -1.2770533585 0.2480275422 0.8950251467 0.2109979248 -1.6249692105 4 0.0826867743 0.7152440948 -1.8916996220 -3.2453517183 -0.2515889959 5 1.3353803719 0.4350300503 0.4779026888 -0.5450539108 0.9029297304 6 1.3265614889 -0.9780238752 1.6815626129 0.9781607829 -0.4352325620 Nov Dec 1 -0.1447615393 0.4989389233 2 0.0265860692 0.0088557808 3 -0.3120808449 -1.7217013050 4 1.2701576671 -0.3167049006 5 -0.1444349622 -0.9627480076 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/16k9g1293630865.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/26k9g1293630865.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/3zc8j1293630865.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/4r3qm1293630865.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/5r3qm1293630865.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/6ndnu1293630865.tab") > > try(system("convert tmp/16k9g1293630865.ps tmp/16k9g1293630865.png",intern=TRUE)) character(0) > try(system("convert tmp/26k9g1293630865.ps tmp/26k9g1293630865.png",intern=TRUE)) character(0) > try(system("convert tmp/3zc8j1293630865.ps tmp/3zc8j1293630865.png",intern=TRUE)) character(0) > try(system("convert tmp/4r3qm1293630865.ps tmp/4r3qm1293630865.png",intern=TRUE)) character(0) > try(system("convert tmp/5r3qm1293630865.ps tmp/5r3qm1293630865.png",intern=TRUE)) character(0) > > > proc.time() user system elapsed 1.473 0.911 4.530