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(1.4,1,-0.8,-2.9,-0.7,-0.7,1.5,3,3.2,3.1,3.9,1,1.3,0.8,1.2,2.9,3.9,4.5,4.5,3.3,2,1.5,1,2.1,3,4,5.1,4.5,4.2,3.3,2.7,1.8,1.4,0.5,-0.4,0.8,0.7,1.9,2,1.1,0.9,0.4,0.7,2.1,2.8,3.9,3.5,2,2,1.5,2.5,3.1,2.7,2.8,2.5,3,3.2,2.8,2.4,2,1.8,1.1,-1.5,-3.7) > 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.27356576 0.47024183 0.00000000 0.08659565 > m$fitted level slope sea Jan 1 1.4000000 0.0000000000 0.000000e+00 Feb 1 1.0380115 -0.1462118136 -1.358882e-02 Mar 1 -0.6847070 -1.1943673594 2.595776e-05 Apr 1 -2.8370746 -1.8215264456 5.980013e-03 May 1 -0.9676461 0.5902089817 -2.404555e-03 Jun 1 -0.6758637 0.3949907020 -2.220614e-03 Jul 1 1.3803564 1.4822045414 -2.431640e-03 Aug 1 2.9928578 1.5674884684 -2.433732e-03 Sep 1 3.2954022 0.7395226926 -2.435418e-03 Oct 1 3.1662779 0.1709519099 -2.436904e-03 Nov 1 3.8637408 0.5155780116 -2.436539e-03 Dec 1 1.2336287 -1.5434227795 -2.437171e-03 Jan 2 0.8511987 -0.8036668065 3.657389e-01 Feb 2 0.7834926 -0.3846271401 -1.808910e-02 Mar 2 1.1647092 0.1201778059 -2.037425e-02 Apr 2 2.8133346 1.1207245948 -2.458387e-02 May 2 3.9252226 1.1149452436 -2.457449e-02 Jun 2 4.5597275 0.8005646371 -2.443111e-02 Jul 2 4.5816065 0.2909113057 -2.438113e-02 Aug 2 3.4303582 -0.6530406549 -2.436814e-02 Sep 2 2.0759175 -1.1121397272 -2.436833e-02 Oct 2 1.4859882 -0.7703284702 -2.436794e-02 Nov 2 1.0032328 -0.5820985108 -2.436785e-02 Dec 2 2.0077592 0.4564213716 -2.436770e-02 Jan 3 2.6843450 0.5986464470 2.996855e-01 Feb 3 3.9788019 1.0150827218 -1.781067e-02 Mar 3 5.1095960 1.0911723404 -1.803835e-02 Apr 3 4.6298944 0.0628897987 -1.517939e-02 May 3 4.2475448 -0.2283861223 -1.486717e-02 Jun 3 3.3629534 -0.6578219397 -1.473784e-02 Jul 3 2.7140805 -0.6519645795 -1.473822e-02 Aug 3 1.8316731 -0.8027992606 -1.473685e-02 Sep 3 1.3883194 -0.5675253480 -1.473679e-02 Oct 3 0.5356905 -0.7541389725 -1.473693e-02 Nov 3 -0.3738424 -0.8558515698 -1.473696e-02 Dec 3 0.6747689 0.3907080158 -1.473684e-02 Jan 4 0.6138485 0.0976578330 1.190561e-01 Feb 4 1.8317279 0.7840970799 4.951663e-04 Mar 4 2.0407001 0.4063491541 1.339306e-03 Apr 4 1.1889728 -0.4171637524 3.049046e-03 May 4 0.8884485 -0.3408479595 2.987982e-03 Jun 4 0.4073007 -0.4326677317 3.008623e-03 Jul 4 0.6475579 0.0077790056 2.987335e-03 Aug 4 1.9983172 0.8868180936 2.981373e-03 Sep 4 2.8030514 0.8330903666 2.981363e-03 Oct 4 3.8791581 0.9921559524 2.981451e-03 Nov 4 3.5911072 0.1542008818 2.981249e-03 Dec 4 2.1167120 -0.9117910009 2.981171e-03 Jan 5 1.9283793 -0.4413236259 1.879506e-02 Feb 5 1.5007348 -0.4328254282 -1.598493e-03 Mar 5 2.4052978 0.4449717775 -3.163778e-03 Apr 5 3.0861495 0.5993726173 -3.419557e-03 May 5 2.7702498 0.0004728159 -3.037266e-03 Jun 5 2.8008285 0.0201763335 -3.040799e-03 Jul 5 2.5248021 -0.1736972650 -3.033324e-03 Aug 5 2.9584025 0.2238058861 -3.035474e-03 Sep 5 3.2016096 0.2365048856 -3.035472e-03 Oct 5 2.8465151 -0.1507245384 -3.035644e-03 Nov 5 2.4230786 -0.3292272218 -3.035678e-03 Dec 5 2.0092532 -0.3846006414 -3.035681e-03 Jan 6 1.7502683 -0.3028050450 4.054741e-02 Feb 6 1.1266504 -0.5039134211 -5.833477e-03 Mar 6 -1.3515127 -1.7991152596 -3.912055e-03 Apr 6 -3.6592815 -2.1320642888 -3.453218e-03 > m$resid Jan Feb Mar Apr May Jun 1 0.000000000 -0.335875850 -1.494946149 -0.901324826 3.510889251 -0.284650174 2 1.232231936 0.560274658 0.725678728 1.449519682 -0.008421064 -0.458427775 3 0.222219739 0.579019444 0.109917033 -1.493032708 -0.424533415 -0.626212395 4 -0.447739861 0.968857689 -0.546995255 -1.197029931 0.111245018 -0.133894783 5 0.710650440 0.012089256 1.272902974 0.224578718 -0.873082653 0.028732514 6 0.122701516 -0.287474333 -1.879963371 -0.484490542 Jul Aug Sep Oct Nov Dec 1 1.585450735 0.124367403 -1.207401631 -0.829132446 0.502559561 -3.002588977 2 -0.743213015 -1.376541329 -0.669492610 0.498454744 0.274491007 1.514447379 3 0.008541629 -0.219958412 0.343094019 -0.272133947 -0.148324918 1.817826437 4 0.642291700 1.281880540 -0.078349790 0.231961336 -1.221968767 -1.554509105 5 -0.282720745 0.579668825 0.018518631 -0.564686913 -0.260305966 -0.080749663 6 > 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/1uzz41259961840.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/27j261259961840.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/3yguz1259961840.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/4808s1259961840.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/5r5z21259961840.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/6udz61259961840.tab") > system("convert tmp/1uzz41259961840.ps tmp/1uzz41259961840.png") > system("convert tmp/27j261259961840.ps tmp/27j261259961840.png") > system("convert tmp/3yguz1259961840.ps tmp/3yguz1259961840.png") > system("convert tmp/4808s1259961840.ps tmp/4808s1259961840.png") > system("convert tmp/5r5z21259961840.ps tmp/5r5z21259961840.png") > > > proc.time() user system elapsed 1.319 0.802 1.595