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(130,136.7,138.1,139.5,140.4,144.6,151.4,147.9,141.5,143.8,143.6,150.5,150.1,154.9,162.1,176.7,186.6,194.8,196.3,228.8,267.2,237.2,254.7,258.2,257.9,269.6,266.9,269.6,253.9,258.6,274.2,301.5,304.5,285.1,287.7,265.5,264.1,276.1,258.9,239.1,250.1,276.8,297.6,295.4,283,275.8,279.7,254.6,234.6,176.9,148.1,122.7,124.9,121.6,128.4,144.5,151.8,167.1,173.8,203.7,199.8) > 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 186.70684 26.75439 0.00000 0.00000 > m$fitted level slope sea Jan 1 130.0000 0.000000000 0.0000000000 Feb 1 136.3671 0.630208185 0.3329421714 Mar 1 137.7687 0.794038098 0.3312754497 Apr 1 139.1697 0.953164950 0.3303079879 May 1 140.0696 0.937815644 0.3303683633 Jun 1 144.2722 1.921949821 0.3277813546 Jul 1 151.0749 3.424130071 0.3251041898 Aug 1 147.5723 1.271171736 0.3277226977 Sep 1 141.1703 -1.124931964 0.3297175775 Oct 1 143.4709 -0.052857323 0.3291057195 Nov 1 143.2709 -0.098962463 0.3291237697 Dec 1 150.1715 2.095107674 0.3285343476 Jan 2 153.3142 2.412692548 -3.2142224955 Feb 2 154.6320 2.090593220 0.2679533168 Mar 2 161.8404 3.699734211 0.2596087239 Apr 2 176.4526 7.125018728 0.2474007209 May 2 186.3547 7.996113278 0.2452684192 Jun 2 194.5548 8.060083832 0.2451609074 Jul 2 196.0525 6.002300254 0.2475350237 Aug 2 228.5590 14.313255878 0.2409531800 Sep 2 266.9632 21.867630197 0.2368466546 Oct 2 236.9571 5.600675405 0.2429161788 Nov 2 254.4580 9.332551332 0.2419604280 Dec 2 257.9577 7.503351665 0.2422819762 Jan 3 261.3060 6.220111848 -3.4060469892 Feb 3 269.2656 6.743858298 0.3343942030 Mar 3 266.5556 3.773882864 0.3444242696 Apr 3 269.2548 3.436654437 0.3452066490 May 3 253.5452 -2.568629193 0.3547735332 Jun 3 258.2477 -0.288398488 0.3522797101 Jul 3 273.8515 4.695163423 0.3485383261 Aug 3 301.1551 11.784894406 0.3448848446 Sep 3 304.1541 9.029713242 0.3458593903 Oct 3 284.7520 0.113526906 0.3480241087 Nov 3 287.3521 0.893333930 0.3478941584 Dec 3 265.1513 -6.349165501 0.3487225694 Jan 4 266.9840 -3.808670815 -2.8840176722 Feb 4 275.7064 0.005714585 0.3935717706 Mar 4 258.4929 -5.401373464 0.4070903329 Apr 4 238.6851 -9.921400101 0.4148520347 May 4 249.6929 -3.357093870 0.4071127214 Jun 4 276.4005 6.071361629 0.3994817075 Jul 4 297.2031 10.690988767 0.3969151917 Aug 4 295.0015 6.647955127 0.3984569783 Sep 4 282.6000 0.674046796 0.4000206589 Oct 4 275.3995 -1.795419124 0.4004643306 Nov 4 279.2998 -0.009226430 0.4002440610 Dec 4 254.1991 -7.878153636 0.4009101137 Jan 5 239.5911 -9.973201625 -4.9910537383 Feb 5 176.8840 -26.126425282 0.0160273089 Mar 5 148.0823 -26.966270680 0.0176940611 Apr 5 122.6830 -26.474683907 0.0170240371 May 5 124.8914 -17.478549288 0.0086058873 Jun 5 121.5943 -13.031149701 0.0057490720 Jul 5 128.3970 -6.811253095 0.0030065604 Aug 5 144.4992 0.374391005 0.0008318168 Sep 5 151.7996 2.546422679 0.0003806057 Oct 5 167.1002 6.546202808 -0.0001897138 Nov 5 173.8002 6.594436508 -0.0001944345 Dec 5 203.7007 13.903485525 -0.0006854294 Jan 6 200.9365 8.706833878 -1.1364513847 > m$resid Jan Feb Mar Apr May Jun 1 0.000000000 0.272842818 0.049983791 0.038081685 -0.003281469 0.199505914 2 0.073332366 -0.055074053 0.309496421 0.660607469 0.168217224 0.012360849 3 -0.271182100 0.094415722 -0.572238443 -0.065092394 -1.160133377 0.440683417 4 0.521294952 0.702245328 -1.042730220 -0.872826881 1.268377484 1.822336847 5 -0.423578759 -3.007739875 -0.162044253 0.094949715 1.738464155 0.859642457 6 -1.041350521 Jul Aug Sep Oct Nov Dec 1 0.296985924 -0.420644306 -0.465548323 0.207751253 -0.008923405 0.424403096 2 -0.397733162 1.606576637 1.460415678 -3.144830212 0.721479913 -0.353639719 3 0.963317664 1.370559629 -0.532643319 -1.723748919 0.150759969 -1.400197394 4 0.893008352 -0.781599823 -1.154912492 -0.477419207 0.345325094 -1.521306646 5 1.202383154 1.389146470 0.419912833 0.773275288 0.009325045 1.413065741 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/1u0271259928671.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/2sflk1259928671.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/39a221259928671.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/4siej1259928671.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/5n5fu1259928671.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/6pvs51259928671.tab") > system("convert tmp/1u0271259928671.ps tmp/1u0271259928671.png") > system("convert tmp/2sflk1259928671.ps tmp/2sflk1259928671.png") > system("convert tmp/39a221259928671.ps tmp/39a221259928671.png") > system("convert tmp/4siej1259928671.ps tmp/4siej1259928671.png") > system("convert tmp/5n5fu1259928671.ps tmp/5n5fu1259928671.png") > > > proc.time() user system elapsed 1.320 0.801 2.222