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(226.9,235.9,216.2,226.2,198.3,176.7,166.2,157.6,163.4,159.7,191.0,239.4,321.9,362.7,413.6,407.1,383.2,347.7,333.8,312.3,295.4,283.3,287.6,265.7,250.2,234.7,244.0,231.2,223.8,223.5,210.5,201.6,190.7,207.5,198.8,196.6,204.2,227.4,229.7,217.9,221.4,216.3,197.0,193.8,196.8,180.5,174.8,181.6,190.0,190.6,179.0,174.1,161.1,168.6,169.4,152.2,148.3,137.7,145.0,153.4) > 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 202.08457 0.00000 22.35167 > m$fitted level slope sea Jan 1 226.9000 0.0000000 0.00000000 Feb 1 234.9776 7.9907545 0.04340841 Mar 1 218.0612 -13.3610226 0.09560943 Apr 1 224.5828 3.8480136 0.07681250 May 1 200.3826 -20.4522262 0.07967290 Jun 1 176.8574 -23.1148359 0.07954858 Jul 1 165.2344 -13.1579002 0.07975667 Aug 1 157.1306 -8.7788576 0.07976649 Sep 1 162.2489 3.2619333 0.07975569 Oct 1 160.0418 -1.4765559 0.07975759 Nov 1 188.6046 24.5501066 0.07975648 Dec 1 237.4476 45.5979369 0.07975684 Jan 2 316.2259 74.0293211 3.14448808 Feb 2 364.7781 54.4036334 -0.61241851 Mar 2 414.5461 50.3657392 -0.59671479 Apr 2 411.6751 4.3237302 -0.50072247 May 2 386.0049 -21.6562516 -0.49351444 Jun 2 349.3507 -34.6516493 -0.49449090 Jul 2 332.8914 -18.8888812 -0.49386149 Aug 2 312.8804 -19.8611504 -0.49386663 Sep 2 295.6881 -17.5487562 -0.49386997 Oct 2 283.3892 -13.0002226 -0.49387345 Nov 2 286.8267 1.2417927 -0.49387489 Dec 2 267.7594 -16.3543990 -0.49387538 Jan 3 245.6630 -21.2935028 4.97645338 Feb 3 234.4272 -13.2936864 -0.35711687 Mar 3 242.7712 5.5280359 -0.41152412 Apr 3 232.7806 -7.8995034 -0.39074066 May 3 224.2400 -8.4548544 -0.39062628 Jun 3 223.3102 -1.9347332 -0.39026260 Jul 3 211.6409 -10.3691054 -0.39051261 Aug 3 201.9391 -9.7909904 -0.39051034 Sep 3 191.1662 -10.6417133 -0.39050943 Oct 3 205.9319 11.3717953 -0.39052192 Nov 3 200.4869 -3.1986521 -0.39052083 Dec 3 197.0118 -3.4381352 -0.39052083 Jan 4 199.1476 1.3652420 4.62506441 Feb 4 225.9594 21.9791779 -0.22925214 Mar 4 231.1713 7.4057779 -0.19573025 Apr 4 219.5370 -9.0733103 -0.17544817 May 4 220.7817 -0.1350742 -0.17691192 Jun 4 216.7755 -3.4892270 -0.17706068 Jul 4 198.3303 -16.4476314 -0.17736611 Aug 4 193.1117 -6.7185351 -0.17733573 Sep 4 196.2198 1.7954668 -0.17734299 Oct 4 181.9182 -12.1513634 -0.17733669 Nov 4 174.6044 -7.9599974 -0.17733694 Dec 4 180.6943 4.2130722 -0.17733674 Jan 5 187.7032 6.6250974 2.08220022 Feb 5 191.0326 3.9234423 -0.20981039 Mar 5 180.3027 -8.8057022 -0.18549701 Apr 5 174.0887 -6.5622128 -0.18778881 May 5 161.7336 -11.5805985 -0.18710673 Jun 5 167.4529 3.4088997 -0.18655495 Jul 5 169.6778 2.3830693 -0.18657502 Aug 5 153.7947 -13.4431707 -0.18661604 Sep 5 147.9044 -6.8992329 -0.18662067 Oct 5 138.1098 -9.4078172 -0.18661973 Nov 5 144.0068 3.8525821 -0.18662038 Dec 5 153.1767 8.4596203 -0.18662032 > m$resid Jan Feb Mar Apr May Jun 1 0.00000000 0.59492982 -1.53209263 1.21454689 -1.70944986 -0.18730274 2 2.12051058 -1.33721830 -0.27870215 -3.22948448 -1.82755058 -0.91415529 3 -0.36124713 0.55126093 1.30551040 -0.94253857 -0.03906599 0.45865576 4 0.34791028 1.42832119 -1.01376894 -1.15724897 0.62875805 -0.23594701 5 0.17368524 -0.18779918 -0.88716734 0.15759542 -0.35301736 1.05443311 Jul Aug Sep Oct Nov Dec 1 0.70042179 0.30804388 0.84700981 -0.33332917 1.83084636 1.48061026 2 1.10883123 -0.06839431 0.16266544 0.31996673 1.00185499 -1.23780463 3 -0.59331570 0.04066752 -0.05984413 1.54854093 -1.02495856 -0.01684645 4 -0.91155878 0.68439357 0.59891772 -0.98109020 0.29484177 0.85631495 5 -0.07216204 -1.11329733 0.46033350 -0.17646644 0.93280319 0.32408224 > 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/1rtcw1260004711.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/2jual1260004711.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/384wa1260004711.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/4o90m1260004711.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/58n5n1260004711.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/6uj211260004711.tab") > system("convert tmp/1rtcw1260004711.ps tmp/1rtcw1260004711.png") > system("convert tmp/2jual1260004711.ps tmp/2jual1260004711.png") > system("convert tmp/384wa1260004711.ps tmp/384wa1260004711.png") > system("convert tmp/4o90m1260004711.ps tmp/4o90m1260004711.png") > system("convert tmp/58n5n1260004711.ps tmp/58n5n1260004711.png") > > > proc.time() user system elapsed 1.261 0.808 2.042