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(5.4,5.4,5.6,5.7,5.8,5.8,5.8,5.9,6.1,6.4,6.4,6.3,6.2,6.2,6.3,6.4,6.5,6.6,6.6,6.6,6.8,7,7.2,7.3,7.5,7.6,7.6,7.7,7.7,7.7,7.7,7.6,7.7,7.9,7.9,7.9,7.8,7.6,7.4,7,7,7.2,7.5,7.8,7.8,7.7,7.6,7.6,7.5,7.5,7.6,7.6,7.9,7.6,7.5,7.5,7.6,7.7,7.8,7.9,7.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.000000e+00 1.760415e-02 5.511737e-05 2.815224e-04 > m$fitted level slope sea Jan 1 5.400000 0.0000000000 0.000000e+00 Feb 1 5.400000 0.0000000000 0.000000e+00 Mar 1 5.594929 0.1867782437 2.214277e-03 Apr 1 5.702526 0.1110168272 -1.395545e-03 May 1 5.800089 0.0981409066 1.029558e-04 Jun 1 5.802352 0.0063765518 -9.821151e-04 Jul 1 5.799926 -0.0020485222 1.999745e-04 Aug 1 5.897437 0.0932387482 1.140485e-03 Sep 1 6.097556 0.1955320400 9.173356e-04 Oct 1 6.397635 0.2955928778 8.713698e-04 Nov 1 6.407470 0.0220980171 -3.387133e-03 Dec 1 6.302380 -0.0996323080 -5.625247e-04 Jan 2 6.199030 -0.1031825136 1.023063e-03 Feb 2 6.200136 -0.0070252437 -1.508804e-03 Mar 2 6.293351 0.0879427324 5.266918e-03 Apr 2 6.401653 0.1071854622 -1.934072e-03 May 2 6.498850 0.0977396072 1.288343e-03 Jun 2 6.599697 0.1006784567 2.599747e-04 Jul 2 6.603385 0.0089552059 -2.042855e-03 Aug 2 6.598814 -0.0038359986 1.372884e-03 Sep 2 6.796094 0.1863574442 1.123048e-03 Oct 2 6.996112 0.1992760596 3.698796e-03 Nov 2 7.201230 0.2048010568 -1.311238e-03 Dec 2 7.302973 0.1073504821 -1.546793e-03 Jan 3 7.497427 0.1896291001 1.368669e-03 Feb 3 7.607032 0.1160160389 -5.989004e-03 Mar 3 7.600149 0.0006306473 1.511380e-03 Apr 3 7.697293 0.0910610195 1.404633e-03 May 3 7.704641 0.0125802444 -3.509246e-03 Jun 3 7.698264 -0.0051943218 1.992635e-03 Jul 3 7.698421 -0.0001774146 1.507122e-03 Aug 3 7.607914 -0.0848654138 -6.692211e-03 Sep 3 7.692333 0.0738470560 5.377450e-03 Oct 3 7.894304 0.1939709893 3.962808e-03 Nov 3 7.905185 0.0223077826 -2.708711e-03 Dec 3 7.905385 0.0015861585 -5.085857e-03 Jan 4 7.801633 -0.0971475181 -2.069863e-04 Feb 4 7.602938 -0.1904762046 -1.620458e-03 Mar 4 7.403631 -0.1987140398 -3.513413e-03 Apr 4 7.004045 -0.3858370403 -1.375219e-03 May 4 6.986309 -0.0427767195 8.787080e-03 Jun 4 7.191292 0.1881634293 5.407114e-03 Jul 4 7.488536 0.2898359474 1.001083e-02 Aug 4 7.809854 0.3191798280 -1.027362e-02 Sep 4 7.816477 0.0278539811 -1.231220e-02 Oct 4 7.694313 -0.1119763663 7.685702e-03 Nov 4 7.601836 -0.0938014118 -2.095377e-03 Dec 4 7.598726 -0.0093051396 6.577127e-05 Jan 5 7.501020 -0.0917198424 1.602843e-04 Feb 5 7.502256 -0.0063154027 -3.458805e-03 Mar 5 7.582862 0.0744182522 1.599008e-02 Apr 5 7.613890 0.0341521150 -1.331854e-02 May 5 7.887817 0.2567254689 9.019253e-03 Jun 5 7.622457 -0.2280094110 -1.556669e-02 Jul 5 7.492263 -0.1371937842 6.445970e-03 Aug 5 7.493302 -0.0088570968 4.874191e-03 Sep 5 7.600379 0.0987785430 -1.909092e-03 Oct 5 7.692374 0.0924800394 7.715773e-03 Nov 5 7.807063 0.1131001730 -7.356508e-03 Dec 5 7.895428 0.0901464749 4.898551e-03 Jan 6 7.907137 0.0172968738 -6.100541e-03 > m$resid Jan Feb Mar Apr May Jun 1 0.00000000 0.00000000 1.42458337 -0.57086225 -0.09704574 -0.69161815 2 -0.02708981 0.72373340 0.71300094 0.14493573 -0.07119348 0.02214978 3 0.62579602 -0.55397450 -0.86786476 0.68095447 -0.59151339 -0.13396466 4 -0.74885179 -0.70241122 -0.06203933 -1.40884140 2.58566112 1.74056265 5 -0.62377921 0.64289291 0.60848594 -0.30314439 1.67752712 -3.65337255 6 -0.55058577 Jul Aug Sep Oct Nov Dec 1 -0.06349889 0.71817008 0.77097372 0.75414795 -2.06130183 -0.91746858 2 -0.69130847 -0.09640595 1.43346786 0.09736624 0.04164138 -0.73448843 3 0.03781190 -0.63828448 1.19619913 0.90536142 -1.29381481 -0.15619258 4 0.76629511 0.22116171 -2.19569223 -1.05388731 0.13698366 0.63698880 5 0.68446762 0.96726000 0.81123846 -0.04747128 0.15541222 -0.17306373 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/1gucg1260533448.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/2m1e31260533448.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/36y421260533448.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/4ykgv1260533448.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/5dy5t1260533448.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/6pzsh1260533448.tab") > system("convert tmp/1gucg1260533448.ps tmp/1gucg1260533448.png") > system("convert tmp/2m1e31260533448.ps tmp/2m1e31260533448.png") > system("convert tmp/36y421260533448.ps tmp/36y421260533448.png") > system("convert tmp/4ykgv1260533448.ps tmp/4ykgv1260533448.png") > system("convert tmp/5dy5t1260533448.ps tmp/5dy5t1260533448.png") > > > proc.time() user system elapsed 1.513 0.824 3.773