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(7.55,7.55,7.59,7.59,7.59,7.57,7.57,7.59,7.6,7.64,7.64,7.76,7.76,7.76,7.77,7.83,7.94,7.94,7.94,8.09,8.18,8.26,8.28,8.28,8.28,8.29,8.3,8.3,8.31,8.33,8.33,8.34,8.48,8.59,8.67,8.67,8.67,8.71,8.72,8.72,8.72,8.74,8.74,8.74,8.74,8.79,8.85,8.86,8.87,8.92,8.96,8.97,8.99,8.98,8.98,9.01,9.01,9.03,9.05,9.05) > 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.0013178681 0.0000417418 0.0000000000 0.0000000000 > m$fitted level slope sea Jan 1 7.550000 0.0000000000 0.000000e+00 Feb 1 7.550000 0.0000000000 0.000000e+00 Mar 1 7.589908 0.0024007106 9.228659e-05 Apr 1 7.589913 0.0021990742 8.721297e-05 May 1 7.589917 0.0019710900 8.304731e-05 Jun 1 7.569954 -0.0006481201 4.638940e-05 Jul 1 7.569953 -0.0005631491 4.732899e-05 Aug 1 7.589927 0.0023154907 7.296657e-05 Sep 1 7.599919 0.0034413709 8.114369e-05 Oct 1 7.639886 0.0089704277 1.141624e-04 Nov 1 7.639893 0.0075833445 1.073133e-04 Dec 1 7.759820 0.0252377562 1.796658e-04 Jan 2 7.768919 0.0228611326 -8.919389e-03 Feb 2 7.759915 0.0177995322 8.548048e-05 Mar 2 7.769912 0.0165450313 8.804051e-05 Apr 2 7.829924 0.0235602888 7.608014e-05 May 2 7.939944 0.0375509980 5.613957e-05 Jun 2 7.939937 0.0314621621 6.339750e-05 Jul 2 7.939932 0.0263541445 6.849128e-05 Aug 2 8.089948 0.0464464389 5.172577e-05 Sep 2 8.179953 0.0535282486 4.678046e-05 Oct 2 8.259956 0.0578344358 4.426366e-05 Nov 2 8.279953 0.0516779852 4.727544e-05 Dec 2 8.279949 0.0432670971 5.071969e-05 Jan 3 8.287818 0.0375929514 -7.818191e-03 Feb 3 8.289807 0.0318899566 1.931701e-04 Mar 3 8.299799 0.0283238665 2.008949e-04 Apr 3 8.299791 0.0237102421 2.092620e-04 May 3 8.309787 0.0214772131 2.126524e-04 Jun 3 8.329787 0.0212366306 2.129582e-04 Jul 3 8.329783 0.0177781396 2.166387e-04 Aug 3 8.339782 0.0165114715 2.177671e-04 Sep 3 8.479797 0.0366210961 2.027687e-04 Oct 3 8.589805 0.0485703759 1.953077e-04 Nov 3 8.669807 0.0536884304 1.926324e-04 Dec 3 8.669804 0.0449457809 1.964582e-04 Jan 4 8.677216 0.0388864246 -7.215580e-03 Feb 4 8.709416 0.0378117009 5.842777e-04 Mar 4 8.719408 0.0332806370 5.918902e-04 Apr 4 8.719400 0.0278593923 5.995154e-04 May 4 8.719395 0.0233217185 6.048587e-04 Jun 4 8.739395 0.0227807245 6.053921e-04 Jul 4 8.739392 0.0190707182 6.084541e-04 Aug 4 8.739389 0.0159650260 6.106000e-04 Sep 4 8.739388 0.0133651659 6.121039e-04 Oct 4 8.789391 0.0193309404 6.092149e-04 Nov 4 8.849393 0.0259535830 6.065300e-04 Dec 4 8.859393 0.0233556878 6.074117e-04 Jan 5 8.877234 0.0224633854 -7.234007e-03 Feb 5 8.919177 0.0256026178 8.231573e-04 Mar 5 8.959180 0.0279480019 8.200111e-04 Apr 5 8.969177 0.0250245636 8.232953e-04 May 5 8.989176 0.0242062132 8.240649e-04 Jun 5 8.979172 0.0186353811 8.284512e-04 Jul 5 8.979170 0.0156005451 8.304517e-04 Aug 5 9.009171 0.0179454774 8.291576e-04 Sep 5 9.009169 0.0150231400 8.305077e-04 Oct 5 9.029170 0.0158335868 8.301943e-04 Nov 5 9.049170 0.0165120512 8.299746e-04 Dec 5 9.049169 0.0138232252 8.307034e-04 > m$resid Jan Feb Mar Apr May Jun 1 0.00000000 0.00000000 1.06827691 -0.06329274 -0.05735050 -0.56800402 2 -0.47956296 -0.68231718 -0.19681318 1.09615309 2.17990158 -0.94682355 3 -0.96993902 -0.81457051 -0.55168991 -0.71385218 -0.34554516 -0.03723108 4 -1.00113850 -0.15726511 -0.70093096 -0.83877536 -0.70215036 -0.08371908 5 -0.14507175 0.46522747 0.36285683 -0.45234836 -0.12663648 -0.86212058 Jul Aug Sep Oct Nov Dec 1 0.01664192 0.52529752 0.19555937 0.92778399 -0.22719616 2.84313202 2 -0.79320012 3.11699661 1.09787823 0.66725949 -0.95364434 -1.30255387 3 -0.53524176 -0.19603854 3.11238489 1.84943285 0.79214922 -1.35315905 4 -0.57415787 -0.48065379 -0.40237966 0.92333941 1.02501873 -0.40209278 5 -0.46968107 0.36292097 -0.45229511 0.12543621 0.10500990 -0.41616853 > 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/1k61q1260535833.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/29lp81260535833.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/3s7l61260535833.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/4f2fv1260535833.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/5jrqy1260535833.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/6ki9x1260535833.tab") > system("convert tmp/1k61q1260535833.ps tmp/1k61q1260535833.png") > system("convert tmp/29lp81260535833.ps tmp/29lp81260535833.png") > system("convert tmp/3s7l61260535833.ps tmp/3s7l61260535833.png") > system("convert tmp/4f2fv1260535833.ps tmp/4f2fv1260535833.png") > system("convert tmp/5jrqy1260535833.ps tmp/5jrqy1260535833.png") > > > proc.time() user system elapsed 1.504 0.857 1.823