R version 2.12.0 (2010-10-15) Copyright (C) 2010 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(103.48,103.93,103.89,104.4,104.79,104.77,105.13,105.26,104.96,104.75,105.01,105.15,105.2,105.77,105.78,106.26,106.13,106.12,106.57,106.44,106.54,107.1,108.1,108.4,108.84,109.62,110.42,110.67,111.66,112.28,112.87,112.18,112.36,112.16,111.49,111.25,111.36,111.74,111.1,111.33,111.25,111.04,110.97,111.31,111.02,111.07,111.36,111.54,112.05,112.52,112.94,113.33,113.78,113.77,113.82,113.89,114.25,114.41) > 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.04710352 0.01879423 0.00000000 0.01433337 > m$fitted level slope sea Jan 1 103.4800 0.000000000 0.00000000 Feb 1 103.8739 0.067154707 0.02051637 Mar 1 103.8787 0.044110707 0.02040372 Apr 1 104.3289 0.212766643 0.01853275 May 1 104.7464 0.299877065 0.01780853 Jun 1 104.7844 0.187752275 0.01834982 Jul 1 105.0963 0.240993750 0.01821155 Aug 1 105.2524 0.204536778 0.01826103 Sep 1 104.9989 0.007804262 0.01839895 Oct 1 104.7622 -0.097258138 0.01843681 Nov 1 104.9553 0.027492066 0.01841375 Dec 1 105.1150 0.084334804 0.01840837 Jan 2 105.3314 0.138908705 -0.14725966 Feb 2 105.7290 0.240710807 0.01946589 Mar 2 105.7839 0.161016128 0.01933369 Apr 2 106.2084 0.274454159 0.01877334 May 2 106.1520 0.132256301 0.01931218 Jun 2 106.1209 0.062124978 0.01946755 Jul 2 106.5098 0.202484849 0.01929992 Aug 2 106.4531 0.091111533 0.01936949 Sep 2 106.5232 0.082118915 0.01937239 Oct 2 107.0278 0.263621837 0.01934228 Nov 2 107.9929 0.565033074 0.01931664 Dec 2 108.4004 0.497342984 0.01931959 Jan 3 109.0146 0.546622951 -0.18890963 Feb 3 109.5986 0.561808998 0.01775163 Mar 3 110.3753 0.654027281 0.01785156 Apr 3 110.6936 0.509582609 0.01831923 May 3 111.5934 0.677266744 0.01790290 Jun 3 112.2630 0.674005487 0.01790763 Jul 3 112.8615 0.641557895 0.01793301 Aug 3 112.3110 0.129413431 0.01814252 Sep 3 112.3528 0.091769414 0.01815047 Oct 3 112.1755 -0.023844525 0.01816303 Nov 3 111.5474 -0.283470880 0.01817750 Dec 3 111.2354 -0.295747606 0.01817785 Jan 4 111.4332 -0.086226374 -0.13421697 Feb 4 111.6882 0.054334513 0.01612804 Mar 4 111.1573 -0.196911244 0.01592582 Apr 4 111.2751 -0.061556660 0.01559988 May 4 111.2321 -0.053588686 0.01558517 Jun 4 111.0415 -0.112453884 0.01564869 Jul 4 110.9515 -0.102784088 0.01564306 Aug 4 111.2448 0.067394000 0.01559131 Sep 4 111.0386 -0.050167142 0.01560978 Oct 4 111.0471 -0.024989508 0.01560774 Nov 4 111.3085 0.098105411 0.01560264 Dec 4 111.5113 0.143072776 0.01560169 Jan 5 112.1181 0.340523531 -0.12555714 Feb 5 112.5033 0.359112168 0.01182444 Mar 5 112.9208 0.384201157 0.01184047 Apr 5 113.3167 0.389217197 0.01183086 May 5 113.7613 0.413003834 0.01179590 Jun 5 113.8043 0.254065492 0.01193239 Jul 5 113.8359 0.158455744 0.01197664 Aug 5 113.8909 0.114040384 0.01198739 Sep 5 114.2121 0.203032499 0.01197627 Oct 5 114.3999 0.196496913 0.01197669 > m$resid Jan Feb Mar Apr May Jun 1 0.00000000 1.05673373 -0.21274522 1.29441439 0.64259368 -0.82006317 2 0.47449731 0.65598954 -0.58260337 0.82273928 -1.03451464 -0.51113546 3 0.39364545 0.10300572 0.67364270 -1.04969404 1.22104152 -0.02377577 4 1.62498585 0.97392946 -1.83463874 0.98458133 0.05804684 -0.42920941 5 1.50853262 0.13030767 0.18316396 0.03650787 0.17333150 -1.15897876 Jul Aug Sep Oct Nov Dec 1 0.38862172 -0.26597607 -1.43510340 -0.76637286 0.90997686 0.41463227 2 1.02359417 -0.81234604 -0.06559437 1.32394345 2.19860238 -0.49375638 3 -0.23664800 -3.73561413 -0.27458631 -0.84332833 -1.89380914 -0.08955096 4 0.07052693 1.24130312 -0.85752731 0.18365457 0.89789932 0.32800854 5 -0.69734880 -0.32397400 0.64913710 -0.04767290 > 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 > postscript(file="/var/www/rcomp/tmp/14jzp1291748448.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') > grid() > dev.off() null device 1 > mylagmax <- nx/2 > postscript(file="/var/www/rcomp/tmp/2fsga1291748448.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/www/rcomp/tmp/3fsga1291748448.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/www/rcomp/tmp/4p1gd1291748448.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/www/rcomp/tmp/5iafy1291748448.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/www/rcomp/createtable file can be downloaded at http://www.wessa.net/cretab > load(file="/var/www/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/rcomp/tmp/6lbd41291748448.tab") > > try(system("convert tmp/14jzp1291748448.ps tmp/14jzp1291748448.png",intern=TRUE)) character(0) > try(system("convert tmp/2fsga1291748448.ps tmp/2fsga1291748448.png",intern=TRUE)) character(0) > try(system("convert tmp/3fsga1291748448.ps tmp/3fsga1291748448.png",intern=TRUE)) character(0) > try(system("convert tmp/4p1gd1291748448.ps tmp/4p1gd1291748448.png",intern=TRUE)) character(0) > try(system("convert tmp/5iafy1291748448.ps tmp/5iafy1291748448.png",intern=TRUE)) character(0) > > > proc.time() user system elapsed 1.430 0.910 2.305