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(9.829,9.125,9.782,9.441,9.162,9.915,10.444,10.209,9.985,9.842,9.429,10.132,9.849,9.172,10.313,9.819,9.955,10.048,10.082,10.541,10.208,10.233,9.439,9.963,10.158,9.225,10.474,9.757,10.490,10.281,10.444,10.640,10.695,10.786,9.832,9.747,10.411,9.511,10.402,9.701,10.540,10.112,10.915,11.183,10.384,10.834,9.886,10.216,10.943,9.867,10.203,10.837,10.573,10.647,11.502,10.656,10.866,10.835,9.945,10.331) > 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') Warning message: In StructTS(x, type = "BSM") : possible convergence problem: optim gave code=1 NEW_X > m$coef level slope seas epsilon 0.0009718043 0.0000000000 0.0361409009 0.0463709171 > m$fitted level slope sea Jan 1 9.829000 0.0000000000 0.000000000 Feb 1 9.808907 -0.0014300831 -0.375425974 Mar 1 9.758719 -0.0116171014 0.111946730 Apr 1 9.665996 -0.0258247874 -0.120154157 May 1 9.523641 -0.0430868496 -0.208772813 Jun 1 9.553772 -0.0337474281 0.255356178 Jul 1 9.718084 -0.0116305272 0.407961734 Aug 1 9.842313 0.0018392640 0.126865418 Sep 1 9.898935 0.0067188303 -0.019005659 Oct 1 9.907491 0.0068673030 -0.069277862 Nov 1 9.839856 0.0013573028 -0.247378593 Dec 1 9.864890 0.0029713283 0.212244730 Jan 2 9.855473 0.0025960122 0.036292296 Feb 2 9.807915 0.0003416808 -0.501015564 Mar 2 9.830876 0.0015309934 0.430675937 Apr 2 9.842602 0.0020819349 -0.045780880 May 2 9.902660 0.0051464435 -0.075840409 Jun 2 9.933583 0.0064517841 0.055223086 Jul 2 9.931226 0.0060280708 0.171878480 Aug 2 9.997643 0.0087790833 0.392704647 Sep 2 10.045591 0.0104672480 0.061049260 Oct 2 10.086623 0.0117106308 0.064595299 Nov 2 10.063544 0.0103884128 -0.528310270 Dec 2 10.027181 0.0087915015 0.070675827 Jan 3 10.028770 0.0086038526 0.151787325 Feb 3 10.009445 0.0078473163 -0.699177293 Mar 3 10.016585 0.0078262287 0.459476113 Apr 3 10.009268 0.0073497508 -0.209250812 May 3 10.074740 0.0092065521 0.251054537 Jun 3 10.117682 0.0102742257 0.067357834 Jul 3 10.163273 0.0113652647 0.178968391 Aug 3 10.198136 0.0120673114 0.373107258 Sep 3 10.261054 0.0135258883 0.282721202 Oct 3 10.325917 0.0149275864 0.304832873 Nov 3 10.350883 0.0151848574 -0.549796930 Dec 3 10.302133 0.0136846890 -0.354020489 Jan 4 10.291523 0.0131772526 0.197794542 Feb 4 10.290470 0.0128833224 -0.733699850 Mar 4 10.267660 0.0121155313 0.247460872 Apr 4 10.245134 0.0113415333 -0.435672124 May 4 10.257201 0.0113580611 0.280539540 Jun 4 10.254140 0.0110289641 -0.097186699 Jul 4 10.312927 0.0121058099 0.452446127 Aug 4 10.392698 0.0135970072 0.576754749 Sep 4 10.401572 0.0134962155 -0.002539435 Oct 4 10.422031 0.0136388136 0.389591735 Nov 4 10.426674 0.0134638746 -0.511460375 Dec 4 10.444381 0.0135414025 -0.242311828 Jan 5 10.480718 0.0139332855 0.386705530 Feb 5 10.504719 0.0141032814 -0.671132978 Mar 5 10.467049 0.0132172044 -0.093047334 Apr 5 10.541170 0.0142812480 0.095733907 May 5 10.550540 0.0141942298 0.038549725 Jun 5 10.591229 0.0146651862 -0.031018374 Jul 5 10.663497 0.0156815490 0.649379979 Aug 5 10.640622 0.0150126071 0.142470325 Sep 5 10.666790 0.0152013284 0.162251743 Oct 5 10.663517 0.0148989297 0.233041862 Nov 5 10.655735 0.0145421398 -0.634693319 Dec 5 10.658632 0.0143666931 -0.288348577 > m$resid Jan Feb Mar Apr May Jun 1 0.00000000 -2.16410084 -0.71112209 -0.90667305 -1.32687523 0.90457491 2 -0.34405496 -1.11272542 0.43644010 0.19046893 1.10284724 0.50802823 3 -0.18870364 -0.71491153 -0.01740237 -0.36496003 1.39637907 0.81639383 4 -0.65946295 -0.38541824 -0.95442976 -0.91706028 0.01913113 -0.38082252 5 0.63708205 0.28161451 -1.44211833 1.68897790 -0.13590308 0.73336305 Jul Aug Sep Oct Nov Dec 1 2.67070561 1.98515729 0.85935703 0.03066679 -1.31335616 0.43782989 2 -0.18044568 1.28309958 0.86015763 0.69163010 -0.81083313 -1.13098957 3 0.86521564 0.58395404 1.28255850 1.31457116 0.26127458 -1.69638232 4 1.26750891 1.80818829 -0.12720081 0.18917564 -0.24672784 0.11753615 5 1.59813697 -1.07373318 0.31211748 -0.51960709 -0.64149111 -0.33120729 > 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/1teo21322255705.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/2nj891322255705.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/36n1l1322255705.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/403j51322255705.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/5undk1322255705.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/6zxl01322255705.tab") > > try(system("convert tmp/1teo21322255705.ps tmp/1teo21322255705.png",intern=TRUE)) character(0) > try(system("convert tmp/2nj891322255705.ps tmp/2nj891322255705.png",intern=TRUE)) character(0) > try(system("convert tmp/36n1l1322255705.ps tmp/36n1l1322255705.png",intern=TRUE)) character(0) > try(system("convert tmp/403j51322255705.ps tmp/403j51322255705.png",intern=TRUE)) character(0) > try(system("convert tmp/5undk1322255705.ps tmp/5undk1322255705.png",intern=TRUE)) character(0) > > > proc.time() user system elapsed 2.228 0.313 2.714