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(8.9,8.8,8.3,7.5,7.2,7.4,8.8,9.3,9.3,8.7,8.2,8.3,8.5,8.6,8.5,8.2,8.1,7.9,8.6,8.7,8.7,8.5,8.4,8.5,8.7,8.7,8.6,8.5,8.3,8,8.2,8.1,8.1,8,7.9,7.9,8,8,7.9,8,7.7,7.2,7.5,7.3,7,7,7,7.2,7.3,7.1,6.8,6.4,6.1,6.5,7.7,7.9,7.5,6.9,6.6,6.9) > 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.0000000000 0.1238038872 0.0009755658 0.0000000000 > m$fitted level slope sea Jan 1 8.900000 0.000000000 0.0000000000 Feb 1 8.801165 -0.098931180 -0.0011651106 Mar 1 8.311165 -0.477189481 -0.0111650143 Apr 1 7.506405 -0.791631310 -0.0064051194 May 1 7.185035 -0.339964679 0.0149648540 Jun 1 7.388638 0.182180592 0.0113617080 Jul 1 8.771030 1.335066581 0.0289704077 Aug 1 9.330217 0.589784159 -0.0302172896 Sep 1 9.311613 0.005386044 -0.0116125968 Oct 1 8.712685 -0.575097458 -0.0126852113 Nov 1 8.194464 -0.520463752 0.0055357217 Dec 1 8.283879 0.065365170 0.0161211387 Jan 2 8.499473 0.209633542 0.0005267953 Feb 2 8.603027 0.107693347 -0.0030272429 Mar 2 8.499393 -0.090630848 0.0006074545 Apr 2 8.221140 -0.266503473 -0.0211400251 May 2 8.068595 -0.159729499 0.0314049335 Jun 2 7.950454 -0.120752516 -0.0504542607 Jul 2 8.523017 0.529022662 0.0769830506 Aug 2 8.741716 0.238192770 -0.0417163437 Sep 2 8.698977 -0.025091816 0.0010227349 Oct 2 8.505567 -0.182837383 -0.0055669077 Nov 2 8.408892 -0.102086749 -0.0088918023 Dec 2 8.484405 0.064354509 0.0155951916 Jan 3 8.694772 0.201158835 0.0052278805 Feb 3 8.707179 0.024173805 -0.0071794132 Mar 3 8.591859 -0.104863339 0.0081410835 Apr 3 8.522675 -0.071813593 -0.0226746433 May 3 8.256898 -0.251282774 0.0431015545 Jun 3 8.097687 -0.166063731 -0.0976874060 Jul 3 8.104620 -0.005937750 0.0953796526 Aug 3 8.130362 0.023384142 -0.0303620022 Sep 3 8.090147 -0.035481564 0.0098529402 Oct 3 8.006559 -0.080008114 -0.0065588684 Nov 3 7.914884 -0.090806451 -0.0148843387 Dec 3 7.896038 -0.024205337 0.0039624488 Jan 4 7.981359 0.077145859 0.0186407720 Feb 4 8.000610 0.023547608 -0.0006096478 Mar 4 7.907995 -0.083175270 -0.0079950276 Apr 4 7.994797 0.073267156 0.0052029328 May 4 7.677933 -0.285370643 0.0220672791 Jun 4 7.312352 -0.359125511 -0.1123516518 Jul 4 7.385142 0.038055289 0.1148578432 Aug 4 7.333887 -0.044071494 -0.0338869911 Sep 4 7.003330 -0.307510243 -0.0033299412 Oct 4 6.985915 -0.040750106 0.0140853458 Nov 4 7.012936 0.021569973 -0.0129355571 Dec 4 7.199209 0.173010844 0.0007908939 Jan 5 7.284697 0.092541883 0.0153031250 Feb 5 7.103319 -0.159316911 -0.0033193556 Mar 5 6.841559 -0.253084869 -0.0415594568 Apr 5 6.372266 -0.451292429 0.0277338183 May 5 6.048116 -0.334852337 0.0518837318 Jun 5 6.605800 0.482676445 -0.1057996861 Jul 5 7.564061 0.918336944 0.1359393346 Aug 5 7.934299 0.416252566 -0.0342993242 Sep 5 7.565976 -0.302447462 -0.0659764403 Oct 5 6.900936 -0.634598887 -0.0009357024 Nov 5 6.607254 -0.322304656 -0.0072540294 Dec 5 6.869279 0.212898058 0.0307208282 > m$resid Jan Feb Mar Apr May Jun 1 0.00000000 -0.28275162 -1.09234441 -0.89287324 1.28373061 1.48396631 2 0.41017163 -0.29088533 -0.56915630 -0.49927505 0.30349446 0.11077449 3 0.38898127 -0.50392109 -0.36858148 0.09386861 -0.51008327 0.24219985 4 0.28820565 -0.15241833 -0.30411173 0.44454001 -1.01918476 -0.20962213 5 -0.22882773 -0.71584364 -0.26688667 -0.56337300 0.33087338 2.32356291 Jul Aug Sep Oct Nov Dec 1 3.27656851 -2.11813549 -1.66089304 -1.64976749 0.15527213 1.66495948 2 1.84670003 -0.82655525 -0.74826994 -0.44832198 0.22949802 0.47303648 3 0.45508717 0.08333452 -0.16729973 -0.12654702 -0.03068956 0.18928488 4 1.12880724 -0.23340919 -0.74870809 0.75814781 0.17711793 0.43040904 5 1.23816733 -1.42695311 -2.04258685 -0.94399437 0.88756063 1.52112584 > 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/1q2421259694227.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/2bj9t1259694227.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/3homh1259694227.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/427u31259694227.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/5w8wc1259694227.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/6q42z1259694227.tab") > > system("convert tmp/1q2421259694227.ps tmp/1q2421259694227.png") > system("convert tmp/2bj9t1259694227.ps tmp/2bj9t1259694227.png") > system("convert tmp/3homh1259694227.ps tmp/3homh1259694227.png") > system("convert tmp/427u31259694227.ps tmp/427u31259694227.png") > system("convert tmp/5w8wc1259694227.ps tmp/5w8wc1259694227.png") > > > proc.time() user system elapsed 1.356 0.829 1.689