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(9487,8700,9627,8947,9283,8829,9947,9628,9318,9605,8640,9214,9567,8547,9185,9470,9123,9278,10170,9434,9655,9429,8739,9552,9687,9019,9672,9206,9069,9788,10312,10105,9863,9656,9295,9946,9701,9049,10190,9706,9765,9893,9994,10433,10073,10112,9266,9820,10097,9115,10411,9678,10408,10153,10368,10581,10597,10680,9738,9556) > 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 5644.591 0.000 48295.254 7758.564 > m$fitted level slope sea Jan 1 9487.000 0.0000000 0.0000000 Feb 1 9410.662 14.8108814 -642.7529153 Mar 1 9348.403 8.0207882 293.1207177 Apr 1 9192.928 -8.5036745 -221.5646229 May 1 9172.063 -9.6570100 113.1561361 Jun 1 9078.494 -16.3160406 -231.2918098 Jul 1 9246.430 -4.2939842 655.5955400 Aug 1 9394.132 3.8878454 194.2867153 Sep 1 9434.334 5.5472260 -126.1184829 Oct 1 9488.819 7.5211549 102.7603094 Nov 1 9323.679 1.1420028 -635.8977447 Dec 1 9234.356 -2.0060517 4.7938659 Jan 2 9230.454 -2.0019639 337.1093875 Feb 2 9186.990 -2.0777977 -627.8441707 Mar 2 9074.752 -4.1549887 139.6934706 Apr 2 9172.045 -0.8865179 273.0174580 May 2 9181.471 -0.4851732 -60.9272221 Jun 2 9308.786 4.6801842 -61.5660557 Jul 2 9429.393 9.1615369 711.8719070 Aug 2 9435.226 9.0443186 -0.3762785 Sep 2 9501.922 10.8372984 138.0190825 Oct 2 9452.069 9.2103888 -6.9508098 Nov 2 9401.891 7.8784491 -646.9312757 Dec 2 9410.676 7.8946876 141.0781807 Jan 3 9395.032 7.5793584 298.4129501 Feb 3 9444.029 8.1194065 -436.3186832 Mar 3 9514.830 9.1257145 140.4158756 Apr 3 9441.238 7.4630639 -213.6474862 May 3 9385.910 5.9944676 -300.7786919 Jun 3 9491.082 8.4899684 271.5283783 Jul 3 9554.950 9.8954018 742.7847143 Aug 3 9698.523 13.1495405 371.6654290 Sep 3 9741.861 13.8276660 113.1871632 Oct 3 9735.314 13.4174295 -73.8942294 Nov 3 9782.620 14.0178772 -496.7022445 Dec 3 9811.360 14.2471906 130.6732111 Jan 4 9751.494 13.2010167 -30.4799534 Feb 4 9684.366 12.0932151 -613.7166589 Mar 4 9732.814 12.6213954 447.4485736 Apr 4 9802.995 13.5334304 -112.2899741 May 4 9921.081 15.3286555 -183.6889484 Jun 4 9929.715 15.2076123 -34.9521925 Jul 4 9820.575 12.9225934 206.2166094 Aug 4 9843.612 13.1051450 586.7095529 Sep 4 9880.143 13.5085258 186.6230459 Oct 4 9964.712 14.6518361 128.2802661 Nov 4 9950.690 14.2251415 -676.9875417 Dec 4 9881.336 13.0701929 -38.8190975 Jan 5 9910.492 13.2808038 182.1685880 Feb 5 9905.377 13.0450385 -785.4170474 Mar 5 9938.632 13.3070737 466.9305017 Apr 5 9935.923 13.0925557 -253.6265830 May 5 10066.419 14.7240774 310.1643861 Jun 5 10125.661 15.3600193 15.4416286 Jul 5 10166.002 15.7206339 195.3191407 Aug 5 10157.492 15.3740445 429.9966973 Sep 5 10218.693 16.0119282 366.0072595 Oct 5 10303.775 16.9350992 357.6400706 Nov 5 10351.257 17.3246717 -621.4952143 Dec 5 10221.099 15.5262153 -625.2561800 > m$resid Jan Feb Mar Apr May Jun 1 0.00000000 -2.62458917 -0.66467133 -1.26282999 -0.11290326 -0.88990626 2 -0.02604033 -0.56230132 -1.37964440 1.19150266 0.11902132 1.49787747 3 -0.30563880 0.53533351 0.79665460 -1.03183614 -0.77488946 1.22357524 4 -0.95529632 -1.03318116 0.46510086 0.73176098 1.32326996 -0.08462745 5 0.20743921 -0.23705906 0.25994883 -0.20553006 1.50389693 0.56987524 Jul Aug Sep Oct Nov Dec 1 2.13703771 1.84896929 0.45265792 0.61755870 -2.19217046 -1.15197052 2 1.39504375 -0.04104841 0.72382383 -0.77123226 -0.76075172 0.01168413 3 0.68799976 1.67702599 0.38225835 -0.25992148 0.43467324 0.18951921 4 -1.57492203 0.12859475 0.29910576 0.91096207 -0.36874396 -1.07702121 5 0.31998447 -0.31085878 0.58910152 0.88970399 0.39417200 -1.90539596 > 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/1bmhq1259939165.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/2l8vb1259939165.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/3hya11259939165.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/4cc6x1259939165.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/5cbwh1259939165.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/6zdck1259939165.tab") > system("convert tmp/1bmhq1259939165.ps tmp/1bmhq1259939165.png") > system("convert tmp/2l8vb1259939165.ps tmp/2l8vb1259939165.png") > system("convert tmp/3hya11259939165.ps tmp/3hya11259939165.png") > system("convert tmp/4cc6x1259939165.ps tmp/4cc6x1259939165.png") > system("convert tmp/5cbwh1259939165.ps tmp/5cbwh1259939165.png") > > > proc.time() user system elapsed 1.608 0.795 1.798