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(115.47,103.34,102.60,100.69,105.67,123.61,113.08,106.46,123.38,109.87,95.74,123.06,123.39,120.28,115.33,110.4,114.49,132.03,123.16,118.82,128.32,112.24,104.53,132.57,122.52,131.8,124.55,120.96,122.6,145.52,118.57,134.25,136.7,121.37,111.63,134.42,137.65,137.86,119.77,130.69,128.28,147.45,128.42,136.9,143.95,135.64,122.48,136.83,153.04,142.71,123.46,144.37,146.15,147.61,158.51,147.4,165.05,154.64,126.2,157.36,154.15,123.21,113.07,110.45,113.57,122.44,114.93,111.85,126.04,121.34) > 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 20.60167 0.00000 16.54628 0.00000 > m$fitted level slope sea Jan 1 115.4700 0.0000000 0.0000000 Feb 1 107.3280 -0.3251796 -3.9879586 Mar 1 103.3357 -0.5311085 -0.7357024 Apr 1 101.4778 -0.5820717 -0.7878003 May 1 103.3585 -0.5305688 2.3115199 Jun 1 114.1760 -0.3868933 9.4339532 Jul 1 115.9997 -0.3637087 -2.9197012 Aug 1 111.6015 -0.4065180 -5.1414575 Sep 1 116.7782 -0.3446806 6.6017558 Oct 1 114.4270 -0.3673630 -4.5569508 Nov 1 104.4747 -0.4756721 -8.7346751 Dec 1 111.9688 -0.3865226 11.0911956 Jan 2 115.7614 -0.4981894 7.6285917 Feb 2 119.2275 -0.5062000 1.0524592 Mar 2 118.0100 -0.5202389 -2.6800400 Apr 2 115.6566 -0.5646565 -5.2566188 May 2 116.2734 -0.5414895 -1.7833810 Jun 2 119.6769 -0.4874003 12.3531026 Jul 2 122.7522 -0.4518220 0.4077730 Aug 2 124.8725 -0.4298863 -6.0524760 Sep 2 122.4705 -0.4461615 5.8495038 Oct 2 117.3938 -0.4824343 -5.1537821 Nov 2 116.4795 -0.4849625 -11.9494933 Dec 2 119.7899 -0.4797834 12.7800756 Jan 3 118.6193 -0.4772479 3.9007217 Feb 3 123.7701 -0.4681738 8.0299498 Mar 3 125.9585 -0.4400294 -1.4084760 Apr 3 127.0911 -0.4155420 -6.1311008 May 3 127.2572 -0.4064605 -4.6571817 Jun 3 131.1269 -0.3504564 14.3931385 Jul 3 125.7807 -0.4022898 -7.2107478 Aug 3 130.7626 -0.3568826 3.4873708 Sep 3 130.2550 -0.3579510 6.4450400 Oct 3 128.3036 -0.3669723 -6.9335958 Nov 3 126.4739 -0.3723450 -14.8438889 Dec 3 124.0587 -0.3754163 10.3612793 Jan 4 129.8661 -0.3703017 7.7839384 Feb 4 131.5376 -0.3632099 6.3224499 Mar 4 127.1652 -0.3949113 -7.3951999 Apr 4 131.5190 -0.3412314 -0.8289770 May 4 133.5613 -0.3117358 -5.2813428 Jun 4 132.7647 -0.3173258 14.6852748 Jul 4 135.1541 -0.2907572 -6.7340768 Aug 4 134.3340 -0.2949985 2.5660075 Sep 4 134.9955 -0.2889260 8.9545216 Oct 4 138.1907 -0.2723977 -2.5506653 Nov 4 138.2545 -0.2713165 -15.7744989 Dec 4 134.1485 -0.2796686 2.6814938 Jan 5 139.0381 -0.2675166 14.0018742 Feb 5 137.9946 -0.2706965 4.7154201 Mar 5 135.6771 -0.2844776 -12.2171027 Apr 5 140.4496 -0.2388820 3.9204484 May 5 146.5271 -0.1750293 -0.3770741 Jun 5 141.4546 -0.2236195 6.1554250 Jul 5 151.9729 -0.1284743 6.5370791 Aug 5 150.6726 -0.1371423 -3.2725784 Sep 5 153.7430 -0.1183526 11.3069882 Oct 5 155.6694 -0.1093170 -1.0294470 Nov 5 148.8381 -0.1312966 -22.6381058 Dec 5 151.7104 -0.1231649 5.6495761 Jan 6 146.3750 -0.1389662 7.7750005 Feb 6 132.1039 -0.1993553 -8.8938796 Mar 6 128.2427 -0.2213964 -15.1726974 Apr 6 118.4373 -0.2943147 -7.9873054 May 6 114.3723 -0.3264131 -0.8022584 Jun 6 117.2689 -0.2988394 5.1710835 Jul 6 113.4846 -0.3262868 1.4454489 Aug 6 114.5146 -0.3171233 -2.6646061 Sep 6 114.7763 -0.3139486 11.2636969 Oct 6 116.8778 -0.3035571 4.4621726 > m$resid Jan Feb Mar Apr May Jun 1 0.00000000 -1.31384092 -0.62699400 -0.27757831 0.54103237 2.50565762 2 0.98764195 0.88447672 -0.14485350 -0.38119865 0.25528299 0.86625354 3 -0.15408501 1.23641176 0.56528748 0.33327916 0.12529895 0.93448053 4 1.36299363 0.44618146 -0.86385854 1.01856100 0.51487377 -0.10577194 5 1.13538098 -0.16943141 -0.44346857 1.09224676 1.36894574 -1.06825724 6 -1.14322065 -3.08655763 -0.79601404 -2.07908279 -0.81952388 0.70345469 Jul Aug Sep Oct Nov Dec 1 0.48778485 -0.88934316 1.22986578 -0.44188343 -2.11065284 1.75506235 2 0.78528059 0.56717229 -0.43477304 -1.02038976 -0.09507371 0.83685262 3 -1.09901390 1.18700408 -0.03324816 -0.35106252 -0.32201849 -0.45012445 4 0.59435383 -0.11656817 0.21074065 0.76720147 0.07399735 -0.84376454 5 2.35575479 -0.25773281 0.70600561 0.44995758 -1.47852021 0.66017909 6 -0.76383774 0.29798288 0.12728244 0.53116410 > 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/1t91j1259862510.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/2j05x1259862510.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/330541259862510.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/4qoel1259862510.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/5kqd51259862510.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/6jp7x1259862510.tab") > > system("convert tmp/1t91j1259862510.ps tmp/1t91j1259862510.png") > system("convert tmp/2j05x1259862510.ps tmp/2j05x1259862510.png") > system("convert tmp/330541259862510.ps tmp/330541259862510.png") > system("convert tmp/4qoel1259862510.ps tmp/4qoel1259862510.png") > system("convert tmp/5kqd51259862510.ps tmp/5kqd51259862510.png") > > > proc.time() user system elapsed 1.372 0.810 1.615