R version 2.13.0 (2011-04-13) Copyright (C) 2011 The R Foundation for Statistical Computing ISBN 3-900051-07-0 Platform: i486-pc-linux-gnu (32-bit) 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.36,1.4,1.4,1.43,1.42,1.44,1.44,1.44,1.42,1.41,1.39,1.41,1.42,1.44,1.43,1.44,1.41,1.43,1.44,1.5,1.59,1.64,1.69,1.8,1.88,1.96,2.04,2.05,2.06,2.04,2.04,2.05,2.01,2.05,2.01,2.03,2.02,2,1.98,1.99,1.99,1.99,1.83,1.76,1.71,1.67,1.65,1.62,1.61,1.6,1.59,1.58,1.58,1.58,1.58,1.58,1.58,1.58,1.58,1.58,1.56,1.62,1.67,1.64,1.63,1.63,1.62) > 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 5.990002e-04 2.283060e-04 7.418454e-07 1.442194e-06 > m$fitted level slope sea Jan 1 1.360000 0.0000000000 0.000000e+00 Feb 1 1.398101 0.0062404446 1.871654e-03 Mar 1 1.398143 0.0037358481 1.865460e-03 Apr 1 1.428119 0.0152440042 1.846438e-03 May 1 1.418269 0.0039665720 1.763400e-03 Jun 1 1.438107 0.0111544121 1.872719e-03 Jul 1 1.438224 0.0061427888 1.790722e-03 Aug 1 1.438186 0.0033338645 1.822294e-03 Sep 1 1.418246 -0.0072469372 1.784363e-03 Oct 1 1.408160 -0.0085376219 1.843303e-03 Nov 1 1.388204 -0.0137298159 1.810507e-03 Dec 1 1.408037 0.0015323767 1.919003e-03 Jan 2 1.432496 0.0116447168 -1.252501e-02 Feb 2 1.438909 0.0094912717 1.095391e-03 Mar 2 1.429169 0.0007078218 8.552208e-04 Apr 2 1.438708 0.0047242217 1.280578e-03 May 2 1.409250 -0.0108039541 7.941509e-04 Jun 2 1.428623 0.0029027191 1.337734e-03 Jul 2 1.438966 0.0062821808 1.024558e-03 Aug 2 1.498576 0.0305073258 1.354941e-03 Sep 2 1.588713 0.0575965193 1.209789e-03 Oct 2 1.638979 0.0542661085 1.030842e-03 Nov 2 1.689043 0.0523569669 9.628667e-04 Dec 2 1.798579 0.0783322158 1.347228e-03 Jan 3 1.889722 0.0840736599 -9.738755e-03 Feb 3 1.959188 0.0778347075 8.262634e-04 Mar 3 2.039646 0.0790289467 3.506128e-04 Apr 3 2.049330 0.0475283950 7.591085e-04 May 3 2.059830 0.0307216293 2.179547e-04 Jun 3 2.039201 0.0074165904 8.656997e-04 Jul 3 2.039365 0.0041253712 6.438749e-04 Aug 3 2.049194 0.0067139014 7.986353e-04 Sep 3 2.009405 -0.0143935685 6.553043e-04 Oct 3 2.048563 0.0099139451 1.367480e-03 Nov 3 2.010074 -0.0120574878 -1.125467e-05 Dec 3 2.028656 0.0018477955 1.304699e-03 Jan 4 2.028935 0.0011419241 -8.933300e-03 Feb 4 1.999783 -0.0120158646 2.494323e-04 Mar 4 1.979415 -0.0158115504 5.958128e-04 Apr 4 1.989108 -0.0042374710 8.595191e-04 May 4 1.989702 -0.0020461788 2.916447e-04 Jun 4 1.989369 -0.0012694903 6.290838e-04 Jul 4 1.831406 -0.0723278331 -1.204693e-03 Aug 4 1.758577 -0.0725550767 1.423315e-03 Sep 4 1.709967 -0.0616957397 2.153587e-06 Oct 4 1.668299 -0.0526124944 1.674784e-03 Nov 4 1.650278 -0.0369240124 -3.222147e-04 Dec 4 1.619011 -0.0343589192 9.817987e-04 Jan 5 1.613304 -0.0214392063 -3.340385e-03 Feb 5 1.599625 -0.0180363952 3.666063e-04 Mar 5 1.589716 -0.0143487023 2.740505e-04 Apr 5 1.579299 -0.0125655992 6.955325e-04 May 5 1.580350 -0.0063954214 -3.671761e-04 Jun 5 1.577877 -0.0046177614 2.117623e-03 Jul 5 1.581166 -0.0010351132 -1.175933e-03 Aug 5 1.578919 -0.0015843847 1.082971e-03 Sep 5 1.580170 -0.0002992016 -1.738425e-04 Oct 5 1.578819 -0.0007761466 1.182799e-03 Nov 5 1.580384 0.0002849827 -3.868354e-04 Dec 5 1.579409 -0.0002857048 5.924991e-04 Jan 6 1.564771 -0.0067630780 -4.752669e-03 Feb 6 1.618797 0.0200556478 1.132092e-03 Mar 6 1.668971 0.0337046181 9.902181e-04 Apr 6 1.639902 0.0052634618 1.781255e-04 May 6 1.629893 -0.0016524836 1.263878e-04 Jun 6 1.627434 -0.0020175416 2.566577e-03 Jul 6 1.620705 -0.0041510782 -6.993180e-04 > m$resid Jan Feb Mar Apr May Jun 1 0.00000000 0.87710652 -0.19329571 0.79660646 -0.75643959 0.47760263 2 0.78400941 -0.12820990 -0.57500135 0.26472382 -1.02651836 0.90683305 3 0.41063782 -0.38991669 0.07855096 -2.07878062 -1.11146450 -1.54203697 4 -0.04913898 -0.83774859 -0.25026649 0.76421608 0.14493963 0.05139434 5 0.88711656 0.21880231 0.24351292 0.11777077 0.40815511 0.11763340 6 -0.44102918 1.73509545 0.90222520 -1.87884545 -0.45751197 -0.02415752 Jul Aug Sep Oct Nov Dec 1 -0.33207150 -0.18596579 -0.70033347 -0.08542300 -0.34363407 1.01008859 2 0.22363798 1.60322661 1.79280740 -0.22041343 -0.12635120 1.71910055 3 -0.21780585 0.17131134 -1.39693123 1.60872266 -1.45411950 0.92029024 4 -4.70256415 -0.01503927 0.71869192 0.60114867 1.03830084 0.16976811 5 0.23709796 -0.03635157 0.08505597 -0.03156526 0.07022815 -0.03777153 6 -0.14119729 > 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/wessaorg/rcomp/tmp/10a7g1322573019.ps",horizontal=F,onefile=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/wessaorg/rcomp/tmp/2d9dt1322573019.ps",horizontal=F,onefile=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/wessaorg/rcomp/tmp/3glr51322573019.ps",horizontal=F,onefile=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/wessaorg/rcomp/tmp/41wrk1322573019.ps",horizontal=F,onefile=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/wessaorg/rcomp/tmp/5b12f1322573019.ps",horizontal=F,onefile=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/wessaorg/rcomp/createtable file can be downloaded at http://www.wessa.net/cretab > load(file="/var/wessaorg/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/wessaorg/rcomp/tmp/6pcmd1322573019.tab") > > try(system("convert tmp/10a7g1322573019.ps tmp/10a7g1322573019.png",intern=TRUE)) character(0) > try(system("convert tmp/2d9dt1322573019.ps tmp/2d9dt1322573019.png",intern=TRUE)) character(0) > try(system("convert tmp/3glr51322573019.ps tmp/3glr51322573019.png",intern=TRUE)) character(0) > try(system("convert tmp/41wrk1322573019.ps tmp/41wrk1322573019.png",intern=TRUE)) character(0) > try(system("convert tmp/5b12f1322573019.ps tmp/5b12f1322573019.png",intern=TRUE)) character(0) > > > proc.time() user system elapsed 1.958 0.245 2.209