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(562,561,555,544,537,543,594,611,613,611,594,595,591,589,584,573,567,569,621,629,628,612,595,597,593,590,580,574,573,573,620,626,620,588,566,557,561,549,532,526,511,499,555,565,542,527,510,514,517,508,493,490,469,478,528,534,518,506,502,516,528) > 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.00000 66.29245 21.45595 0.00000 > m$fitted level slope sea Jan 1 562.0000 0.0000000 0.0000000 Feb 1 561.3396 -0.6884348 -0.3396342 Mar 1 556.3017 -3.5375074 -1.3016700 Apr 1 546.1728 -8.0927262 -2.1727609 May 1 536.7357 -9.0118004 0.2642758 Jun 1 539.0494 -1.3776720 3.9506149 Jul 1 580.4943 27.5883281 13.5056914 Aug 1 613.8826 31.5182884 -2.8826298 Sep 1 621.9885 15.6535592 -8.9884781 Oct 1 616.1152 1.0696919 -5.1151595 Nov 1 598.3750 -11.6711827 -4.3749655 Dec 1 591.4653 -8.4460837 3.5346889 Jan 2 589.0139 -4.3987711 1.9861094 Feb 2 587.2956 -2.5812119 1.7043719 Mar 2 582.4630 -4.0626675 1.5370452 Apr 2 573.1471 -7.4818028 -0.1471425 May 2 568.3952 -5.6847704 -1.3952440 Jun 2 576.3040 3.1972589 -7.3039656 Jul 2 604.1093 19.2226737 16.8906786 Aug 2 626.7403 21.4465947 2.2596567 Sep 2 635.0493 12.8657686 -7.0493028 Oct 2 619.8249 -5.4820767 -7.8249425 Nov 2 602.8700 -12.9721762 -7.8700323 Dec 2 593.5588 -10.5849500 3.4412148 Jan 3 589.6470 -6.2324827 3.3530000 Feb 3 585.7351 -4.7169355 4.2649472 Mar 3 576.1282 -7.8825431 3.8718129 Apr 3 571.6796 -5.6705277 2.3203902 May 3 576.3876 1.0503695 -3.3875729 Jun 3 588.0215 7.8982255 -15.0214806 Jul 3 604.2966 13.3000683 15.7034455 Aug 3 619.4823 14.5163407 6.5176639 Sep 3 621.4006 6.3814981 -1.4005710 Oct 3 600.4502 -11.2715551 -12.4502341 Nov 3 577.7505 -18.6492876 -11.7505202 Dec 3 556.4467 -20.3621949 0.5533228 Jan 4 551.7450 -10.2490924 9.2550354 Feb 4 542.1758 -9.8101513 6.8242478 Mar 4 528.9839 -11.9838494 3.0160613 Apr 4 524.1042 -7.4280374 1.8958191 May 4 518.8407 -6.0362797 -7.8407308 Jun 4 517.6543 -2.9159789 -18.6543130 Jul 4 535.2698 10.2673767 19.7302446 Aug 4 553.0046 15.0590465 11.9953797 Sep 4 542.3899 -1.4245081 -0.3898525 Oct 4 534.9573 -5.2830472 -7.9573125 Nov 4 523.4861 -9.2565906 -13.4860833 Dec 4 517.5875 -7.1002426 -3.5874588 Jan 5 508.2411 -8.5434322 8.7588714 Feb 5 499.1427 -8.8998402 8.8573476 Mar 5 490.8340 -8.5212370 2.1660192 Apr 5 486.3285 -5.9526768 3.6714976 May 5 480.0726 -6.1469146 -11.0726044 Jun 5 497.7028 9.0970186 -19.7028018 Jul 5 511.0389 11.8118521 16.9610727 Aug 5 515.5051 7.1118917 18.4948559 Sep 5 517.4111 3.7809184 0.5889362 Oct 5 514.3697 -0.5851492 -8.3696842 Nov 5 515.8006 0.7051952 -13.8005976 Dec 5 518.2319 1.8102893 -2.2318844 Jan 6 518.8688 1.0588277 9.1311607 > m$resid Jan Feb Mar Apr May Jun 1 0.00000000 -0.10293939 -0.36179835 -0.59013091 -0.11173249 0.93718499 2 0.49887838 0.22685583 -0.18050232 -0.42644145 0.22165648 1.08661120 3 0.53627317 0.18644341 -0.38766921 0.27285920 0.82879434 0.83967547 4 1.24420720 0.05390675 -0.26660058 0.56057722 0.17137717 0.38310900 5 -0.17739663 -0.04375566 0.04646258 0.31579632 -0.02389764 1.87259510 6 -0.09232752 Jul Aug Sep Oct Nov Dec 1 3.56271197 0.48265341 -1.94848619 -1.79122630 -1.56482869 0.39610553 2 1.96869232 0.27331676 -1.05381570 -2.25358415 -0.91996412 0.29332941 3 0.66276193 0.14942955 -0.99920324 -2.16813124 -0.90619476 -0.21053890 4 1.61746664 0.58841123 -2.02469279 -0.47393614 -0.48814126 0.26503569 5 0.33317885 -0.57701095 -0.40911387 -0.53633741 0.15854158 0.13581260 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/1its81259934294.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/2agfo1259934294.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/3y5l51259934294.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/4kada1259934294.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/5qkyx1259934294.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/6r75t1259934294.tab") > system("convert tmp/1its81259934294.ps tmp/1its81259934294.png") > system("convert tmp/2agfo1259934294.ps tmp/2agfo1259934294.png") > system("convert tmp/3y5l51259934294.ps tmp/3y5l51259934294.png") > system("convert tmp/4kada1259934294.ps tmp/4kada1259934294.png") > system("convert tmp/5qkyx1259934294.ps tmp/5qkyx1259934294.png") > > > proc.time() user system elapsed 1.254 0.805 1.461