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(6.3,6.2,6.1,6.3,6.5,6.6,6.5,6.2,6.2,5.9,6.1,6.1,6.1,6.1,6.1,6.4,6.7,6.9,7,7,6.8,6.4,5.9,5.5,5.5,5.6,5.8,5.9,6.1,6.1,6,6,5.9,5.5,5.6,5.4,5.2,5.2,5.2,5.5,5.8,5.8,5.5,5.3,5.1,5.2,5.8,5.8,5.5,5,4.9,5.3,6.1,6.5,6.8,6.6,6.4,6.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.0000000000 0.0597732088 0.0001192175 0.0000000000 > m$fitted level slope sea Jan 1 6.300000 0.0000000000 0.000000e+00 Feb 1 6.200297 -0.0997273914 -2.971680e-04 Mar 1 6.099795 -0.1004952434 2.048653e-04 Apr 1 6.297754 0.1945004423 2.245729e-03 May 1 6.500562 0.2027136655 -5.624460e-04 Jun 1 6.600852 0.1014539069 -8.524269e-04 Jul 1 6.501395 -0.0971738506 -1.394867e-03 Aug 1 6.201220 -0.2978670604 -1.220378e-03 Sep 1 6.197397 -0.0071651878 2.602544e-03 Oct 1 5.902826 -0.2913042499 -2.826393e-03 Nov 1 6.095773 0.1874421292 4.227219e-03 Dec 1 6.102376 0.0086583262 -2.375588e-03 Jan 2 6.099641 -0.0026047007 3.594160e-04 Feb 2 6.098098 -0.0015542059 1.901688e-03 Mar 2 6.102915 0.0046948671 -2.915346e-03 Apr 2 6.394245 0.2852705138 5.755009e-03 May 2 6.701308 0.3066062476 -1.308362e-03 Jun 2 6.902127 0.2030341333 -2.126502e-03 Jul 2 6.997768 0.0978920113 2.231668e-03 Aug 2 7.008316 0.0123777494 -8.315922e-03 Sep 2 6.793266 -0.2102843931 6.733782e-03 Oct 2 6.411343 -0.3783269717 -1.134286e-02 Nov 2 5.895190 -0.5132650427 4.809967e-03 Dec 2 5.500195 -0.3974737287 -1.945163e-04 Jan 3 5.494657 -0.0137777566 5.342806e-03 Feb 3 5.593731 0.0967271380 6.269097e-03 Mar 3 5.806747 0.2098520749 -6.746684e-03 Apr 3 5.898523 0.0951171887 1.476710e-03 May 3 6.097779 0.1963158443 2.221152e-03 Jun 3 6.104331 0.0119001182 -4.330513e-03 Jul 3 5.999914 -0.1011371117 8.553866e-05 Aug 3 6.002964 0.0001126418 -2.964331e-03 Sep 3 5.894346 -0.1055525779 5.653525e-03 Oct 3 5.512696 -0.3738667866 -1.269608e-02 Nov 3 5.583337 0.0581108706 1.666257e-02 Dec 3 5.413024 -0.1638733680 -1.302376e-02 Jan 4 5.195903 -0.2156152324 4.096649e-03 Feb 4 5.194281 -0.0076199543 5.718845e-03 Mar 4 5.201775 0.0069965707 -1.775041e-03 Apr 4 5.498503 0.2870710012 1.496500e-03 May 4 5.795352 0.2965218904 4.648185e-03 Jun 4 5.805395 0.0195907459 -5.395053e-03 Jul 4 5.510502 -0.2844114499 -1.050199e-02 Aug 4 5.302001 -0.2110310540 -2.001296e-03 Sep 4 5.078877 -0.2227208818 2.112268e-02 Oct 4 5.221263 0.1302159020 -2.126324e-02 Nov 4 5.770430 0.5352033869 2.956984e-02 Dec 4 5.820970 0.0666979291 -2.096960e-02 Jan 5 5.508994 -0.2993241335 -8.993502e-03 Feb 5 4.991928 -0.5098345194 8.071505e-03 Mar 5 4.898901 -0.1083733119 1.099191e-03 Apr 5 5.292973 0.3755653390 7.026778e-03 May 5 6.084119 0.7757757081 1.588069e-02 Jun 5 6.504602 0.4335983372 -4.601920e-03 Jul 5 6.816739 0.3166204968 -1.673880e-02 Aug 5 6.605884 -0.1913807604 -5.884199e-03 Sep 5 6.380683 -0.2239525719 1.931701e-02 Oct 5 6.435975 0.0449831726 -3.597540e-02 > m$resid Jan Feb Mar Apr May Jun 1 0.000000000 -0.408489296 -0.003154559 1.206502729 0.033593944 -0.414174751 2 -0.046073138 0.004302052 0.025657839 1.147391406 0.087268661 -0.423632791 3 1.569578202 0.452416956 0.464098493 -0.469164425 0.413933193 -0.754300675 4 -0.211663424 0.851278106 0.059918714 1.145248536 0.038657232 -1.132708458 5 -1.497342062 -0.861346784 1.644738648 1.978940984 1.636986656 -1.399580134 Jul Aug Sep Oct Nov Dec 1 -0.812431313 -0.820879469 1.189034741 -1.162191400 1.958178225 -0.731265165 2 -0.430054458 -0.349772182 -0.910737246 -0.687331190 -0.551926460 0.473612060 3 -0.462347214 0.414133810 -0.432194041 -1.097464263 1.766883949 -0.907964943 4 -1.243436196 0.300141992 -0.047813973 1.443589252 1.656488482 -1.916292087 5 -0.478465254 -2.077837110 -0.133225891 1.100006480 > 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/1ydlx1259941783.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/2za7y1259941783.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/35roq1259941783.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/4ezht1259941783.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/5miob1259941783.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/6zxbz1259941783.tab") > system("convert tmp/1ydlx1259941783.ps tmp/1ydlx1259941783.png") > system("convert tmp/2za7y1259941783.ps tmp/2za7y1259941783.png") > system("convert tmp/35roq1259941783.ps tmp/35roq1259941783.png") > system("convert tmp/4ezht1259941783.ps tmp/4ezht1259941783.png") > system("convert tmp/5miob1259941783.ps tmp/5miob1259941783.png") > > > proc.time() user system elapsed 1.351 0.801 1.956