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(441,449,452,462,455,461,461,463,462,456,455,456,472,472,471,465,459,465,468,467,463,460,462,461,476,476,471,453,443,442,444,438,427,424,416,406,431,434,418,412,404,409,412,406,398,397,385,390,413,413,401,397,397,409,419,424,428,430,424,433,456) > 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 9.189080 1.978337 6.409729 0.000000 > m$fitted level slope sea Jan 1 441.0000 0.00000000 0.00000000 Feb 1 446.7174 0.63158152 2.28256424 Mar 1 450.9369 1.33091070 1.06311511 Apr 1 458.5059 3.01370916 3.49405818 May 1 457.9759 1.94881762 -2.97586834 Jun 1 459.8257 1.91822347 1.17434694 Jul 1 460.9487 1.67009185 0.05132375 Aug 1 462.4535 1.61824530 0.54648919 Sep 1 462.4158 1.09724761 -0.41581048 Oct 1 458.1696 -0.58604205 -2.16962782 Nov 1 454.9120 -1.42807380 0.08802017 Dec 1 454.4357 -1.12802594 1.56428966 Jan 2 466.5397 2.93352447 5.46033929 Feb 2 472.3331 3.82385844 -0.33306805 Mar 2 474.1691 3.22512895 -3.16913170 Apr 2 466.9192 0.06277189 -1.91921363 May 2 463.2634 -1.07862675 -4.26340158 Jun 2 463.1298 -0.78744558 1.87020844 Jul 2 466.3048 0.43084587 1.69523434 Aug 2 466.6378 0.40079574 0.36221780 Sep 2 463.3949 -0.71904181 -0.39491317 Oct 2 461.3356 -1.13112980 -1.33558932 Nov 2 462.1296 -0.53990698 -0.12961890 Dec 2 464.0769 0.22184823 -3.07688940 Jan 3 468.8572 1.61929019 7.14278967 Feb 3 473.1858 2.45118923 2.81417339 Mar 3 471.8762 1.31016887 -0.87618374 Apr 3 460.6324 -2.48167323 -7.63241060 May 3 450.9488 -4.67146759 -7.94882219 Jun 3 442.9679 -5.68241377 -0.96788487 Jul 3 439.8035 -4.91347898 4.19648938 Aug 3 435.6622 -4.67800505 2.33780068 Sep 3 428.4721 -5.44363399 -1.47207545 Oct 3 424.9517 -4.85774232 -0.95171307 Nov 3 419.0492 -5.17571485 -3.04918045 Dec 3 412.5036 -5.59254500 -6.50357345 Jan 4 417.8721 -2.25051121 13.12793666 Feb 4 424.5186 0.46198140 9.48137965 Mar 4 418.1151 -1.61935057 -0.11510288 Apr 4 415.7358 -1.84891811 -3.73583839 May 4 412.0749 -2.39769755 -8.07493888 Jun 4 410.4709 -2.15649410 -1.47092311 Jul 4 408.0781 -2.22836911 3.92192898 Aug 4 403.7262 -2.87334762 2.27379267 Sep 4 399.9151 -3.15780733 -1.91511524 Oct 4 396.8745 -3.12230059 0.12553168 Nov 4 390.6210 -4.07085663 -5.62102744 Dec 4 395.9619 -1.21733590 -5.96188641 Jan 5 400.7251 0.59875777 12.27487350 Feb 5 401.3228 0.59842229 11.67722094 Mar 5 401.1358 0.36080833 -0.13582101 Apr 5 400.6343 0.10053472 -3.63428833 May 5 403.6285 0.97482063 -6.62854137 Jun 5 408.2943 2.09257716 0.70573991 Jul 5 412.8736 2.84632530 6.12641155 Aug 5 419.1845 3.89555508 4.81547925 Sep 5 427.3415 5.18422601 0.65852424 Oct 5 431.0440 4.73662142 -1.04400784 Nov 5 434.7107 4.41344429 -10.71074257 Dec 5 440.4686 4.81997614 -7.46864588 Jan 6 444.2400 4.50252967 11.76000385 > m$resid Jan Feb Mar Apr May 1 0.0000000000 1.2855203256 0.8247865203 1.5055713436 -0.8100411953 2 3.1142989379 0.6287547872 -0.4126711180 -2.2785760209 -0.8188780033 3 1.0114774404 0.5912858474 -0.8009514330 -2.7029558090 -1.5697656468 4 2.3923850721 1.9259584805 -1.4702326819 -0.1632361959 -0.3923071557 5 1.2954848725 -0.0002381533 -0.1682668253 -0.1849766590 0.6238367802 6 -0.2261208886 Jun Jul Aug Sep Oct 1 -0.0221539072 -0.1777451310 -0.0370003040 -0.3710618366 -1.1976601523 2 0.2060528347 0.8611522851 -0.0213084369 -0.7955441082 -0.2929048315 3 -0.7196304902 0.5448742107 0.1669060828 -0.5434626630 0.4163307341 4 0.1720112857 -0.0510428035 -0.4573338118 -0.2018229939 0.0252309016 5 0.7971579927 0.5358992698 0.7444152529 0.9143116662 -0.3180869436 6 Nov Dec 1 -0.5988389647 0.2133498281 2 0.4203997317 0.5441637494 3 -0.2263521482 -0.2978846652 4 -0.6756193998 2.0371767761 5 -0.2301846513 0.2899502252 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/1iv0z1261004281.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/257gz1261004281.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/3w62o1261004281.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/4y3t81261004281.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/5zefx1261004281.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/62gxi1261004281.tab") > try(system("convert tmp/1iv0z1261004281.ps tmp/1iv0z1261004281.png",intern=TRUE)) character(0) > try(system("convert tmp/257gz1261004281.ps tmp/257gz1261004281.png",intern=TRUE)) character(0) > try(system("convert tmp/3w62o1261004281.ps tmp/3w62o1261004281.png",intern=TRUE)) character(0) > try(system("convert tmp/4y3t81261004281.ps tmp/4y3t81261004281.png",intern=TRUE)) character(0) > try(system("convert tmp/5zefx1261004281.ps tmp/5zefx1261004281.png",intern=TRUE)) character(0) > > > proc.time() user system elapsed 1.313 0.788 2.112