R version 2.8.0 (2008-10-20) Copyright (C) 2008 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. Natural language support but running in an English locale 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(464,460,467,460,448,443,436,431,484,510,513,503,471,471,476,475,470,461,455,456,517,525,523,519,509,512,519,517,510,509,501,507,569,580,578,565,547,555,562,561,555,544,537,543,594,611,613,611,594,595,591,589,584,573,567,569,621,629,628,612,595,597,593,590,580,574,573,573,620,626,620,588,566,557,561,549,532,526,511,499,555,565,542,527,510,514,517,508,493,490,469,478,528,534,518,506,502,516,528,533,536,537,524,536,587,597,581,564,558,575,580,575,563,552,537,545,601,604,586,564,549) > 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 2.282846 60.760380 11.936393 0.000000 > m$fitted level slope sea Jan 1 464.0000 0.00000000 0.0000000 Feb 1 460.8725 -2.89547459 -0.8724928 Mar 1 464.5340 1.66348640 2.4659898 Apr 1 461.6877 -1.62908510 -1.6877477 May 1 450.0883 -8.77442844 -2.0882673 Jun 1 441.9253 -8.33706936 1.0746569 Jul 1 435.4367 -7.01129814 0.5632759 Aug 1 430.4806 -5.53666978 0.5193636 Sep 1 471.7032 28.00733981 12.2968269 Oct 1 510.7085 35.89576962 -0.7085197 Nov 1 520.6001 17.24466511 -7.6001270 Dec 1 508.5166 -3.79106132 -5.5166036 Jan 2 476.9558 -23.61004013 -5.9558162 Feb 2 468.5748 -12.81115860 2.4252421 Mar 2 470.3696 -2.80071813 5.6304344 Apr 2 473.4917 1.29398590 1.5082727 May 2 472.4864 -0.30004779 -2.4863989 Jun 2 460.2766 -8.50552497 0.7233555 Jul 2 451.5042 -8.68960679 3.4958183 Aug 2 463.0130 5.25919542 -7.0129835 Sep 2 502.0167 28.56858101 14.9832857 Oct 2 524.7015 24.50465742 0.2985010 Nov 2 528.3335 10.09640135 -5.3335209 Dec 2 519.8409 -2.71864903 -0.8409496 Jan 3 517.7368 -2.29467521 -8.7367515 Feb 3 514.8243 -2.72001840 -2.8242678 Mar 3 515.2499 -0.57943034 3.7500962 Apr 3 514.6341 -0.60423492 2.3659113 May 3 508.2325 -4.57314718 1.7674969 Jun 3 505.0360 -3.63381599 3.9640334 Jul 3 503.1771 -2.42386162 -2.1771142 Aug 3 520.3748 10.96872855 -13.3747957 Sep 3 549.5845 23.42662672 19.4155311 Oct 3 574.2542 24.27551733 5.7457791 Nov 3 582.3401 13.22873834 -4.3400731 Dec 3 571.9386 -2.88616962 -6.9386435 Jan 4 558.0937 -10.36707014 -11.0936917 Feb 4 556.5593 -4.34959108 -1.5592790 Mar 4 557.1554 -0.99921273 4.8445599 Apr 4 555.9782 -1.11988512 5.0218071 May 4 552.2211 -2.91284824 2.7789171 Jun 4 541.1763 -8.43431173 2.8236700 Jul 4 542.4263 -1.86763504 -5.4262766 Aug 4 559.0182 10.65711493 -16.0181893 Sep 4 575.7508 14.78202739 18.2491935 Oct 4 598.7372 20.35185031 12.2627945 Nov 4 611.9423 15.50263951 1.0576933 Dec 4 617.0223 8.43062704 -6.0222816 Jan 5 611.6451 -0.94397577 -17.6450896 Feb 5 600.8296 -7.63495682 -5.8295797 Mar 5 587.4147 -11.53976930 3.5852542 Apr 5 580.4564 -8.44387209 8.5435961 May 5 576.1497 -5.64234888 7.8503471 Jun 5 571.9141 -4.69007526 1.0859297 Jul 5 575.4617 0.87975421 -8.4616767 Aug 5 585.4552 7.04241106 -16.4552367 Sep 5 603.7149 14.63089682 17.2851247 Oct 5 616.6692 13.49672001 12.3307530 Nov 5 625.7163 10.48719192 2.2837006 Dec 5 618.7659 -1.30959144 -6.7658682 Jan 6 610.6499 -5.91554618 -15.6499269 Feb 6 601.5512 -8.06672537 -4.5511922 Mar 6 590.9408 -9.78206121 2.0591565 Apr 6 581.4800 -9.56541593 8.5200042 May 6 571.5682 -9.79940763 8.4318339 Jun 6 571.9553 -2.91843837 2.0447133 Jul 6 581.5783 5.54561390 -8.5783007 Aug 6 592.6620 9.28215094 -19.6619651 Sep 6 603.2068 10.13427932 16.7931828 Oct 6 613.0524 9.93943574 12.9476076 Nov 6 614.9519 4.51325544 5.0481423 Dec 6 598.3732 -9.72624927 -10.3732479 Jan 7 581.7816 -14.36164564 -15.7815845 Feb 7 561.7906 -18.15818390 -4.7905830 Mar 7 554.7628 -10.66270741 6.2371985 Apr 7 540.6180 -13.00787738 8.3820327 May 7 526.2034 -13.95630871 5.7965698 Jun 7 524.4163 -5.74957568 1.5836518 Jul 7 521.4851 -3.85027058 -10.4850948 Aug 7 519.6861 -2.46855325 -20.6860611 Sep 7 534.8952 9.43911160 20.1048160 Oct 7 549.4142 12.86128750 15.5857693 Nov 7 537.2153 -4.02399940 4.7847110 Dec 7 533.7141 -3.67158598 -6.7141465 Jan 8 525.7663 -6.55374828 -15.7662976 Feb 8 521.3148 -5.13821453 -7.3147979 Mar 8 510.4404 -8.99637104 6.5595577 Apr 8 499.0081 -10.63497749 8.9919493 May 8 489.7826 -9.68604000 3.2174385 Jun 8 486.7163 -5.22804011 3.2836687 Jul 8 480.9975 -5.55834842 -11.9975381 Aug 8 499.1504 10.39194510 -21.1503628 Sep 8 510.2849 10.89138091 17.7151289 Oct 8 512.8374 5.28216759 21.1626401 Nov 8 514.2502 2.67871508 3.7497721 Dec 8 512.9261 -0.01545987 -6.9261160 Jan 9 517.2954 2.93544879 -15.2954175 Feb 9 521.1730 3.56897910 -5.1730194 Mar 9 520.9939 1.05084493 7.0060511 Apr 9 523.1977 1.82550430 9.8022650 May 9 531.9128 6.45791241 4.0872133 Jun 9 535.0382 4.21655393 1.9617897 Jul 9 542.8578 6.63882554 -18.8578280 Aug 9 556.0767 11.05983101 -20.0766535 Sep 9 567.0125 10.97656104 19.9874823 Oct 9 575.1030 9.03771344 21.8970405 Nov 9 577.6812 4.69683905 3.3188420 Dec 9 574.3610 -0.69234686 -10.3609957 Jan 10 573.2540 -0.97105596 -15.2539736 Feb 10 577.4221 2.48051460 -2.4221304 Mar 10 574.8549 -0.90724825 5.1450871 Apr 10 568.1473 -4.80040647 6.8527265 May 10 558.6044 -7.98543155 4.3956003 Jun 10 551.1231 -7.64673959 0.8768690 Jul 10 555.3175 0.30513779 -18.3175481 Aug 10 563.8647 5.83678990 -18.8647344 Sep 10 578.3273 11.62417932 22.6726719 Oct 10 582.2193 6.43584504 21.7806659 Nov 10 581.6148 1.71000104 4.3851772 Dec 10 576.4525 -2.90439234 -12.4524851 Jan 11 567.7126 -6.82188813 -18.7126182 > m$resid Jan Feb Mar Apr May 1 0.000000000 -0.437932872 0.615780668 -0.427947674 -0.909693445 2 -2.577257209 1.392884431 1.284693579 0.531794164 -0.203882936 3 0.054637584 -0.054529743 0.274414714 -0.003200666 -0.509255256 4 -0.961637941 0.771212589 0.429606186 -0.015528384 -0.230222161 5 -1.203745802 -0.857654463 -0.500792091 0.397897530 0.359741723 6 -0.591113508 -0.275781652 -0.220015703 0.027827142 -0.030044078 7 -0.594729720 -0.486772820 0.961460327 -0.301126612 -0.121764744 8 -0.369733101 0.181507543 -0.494912773 -0.210360455 0.121819365 9 0.378522908 0.081239702 -0.323026516 0.099436524 0.594641092 10 -0.035749538 0.442627136 -0.434589382 -0.499689453 -0.408822994 11 -0.502481685 Jun Jul Aug Sep Oct 1 0.056166256 0.170125038 0.189170319 4.303392776 1.012000308 2 -1.050749329 -0.023629796 1.789594149 2.990216998 -0.521391512 3 0.120258415 0.155220815 1.718473666 1.598166426 0.108906741 4 -0.707371001 0.842033576 1.606965560 0.529188180 0.714553918 5 0.122070389 0.714112952 0.790564126 0.973554356 -0.145511997 6 0.882396407 1.085210408 0.479274286 0.109322827 -0.025000063 7 1.052650488 0.243535428 0.177214938 1.527666675 0.439126585 8 0.571890040 -0.042356487 2.045654311 0.064073230 -0.719799275 9 -0.287552699 0.310635771 0.566993203 -0.010682691 -0.248808417 10 0.043453986 1.019814109 0.709433421 0.742456891 -0.665814771 11 Nov Dec 1 -2.392732469 -2.698654877 2 -1.848420754 -1.645759780 3 -1.417336332 -2.069991864 4 -0.622255654 -0.908207096 5 -0.386221867 -1.514575326 6 -0.696395793 -1.827815357 7 -2.167117637 0.045229989 8 -0.334139738 -0.345745614 9 -0.557127537 -0.691552085 10 -0.606529312 -0.592101856 11 > 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/freestat/rcomp/tmp/1fl671293483851.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/html/freestat/rcomp/tmp/2qc5t1293483851.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/html/freestat/rcomp/tmp/31mmv1293483851.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/html/freestat/rcomp/tmp/41mmv1293483851.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/www/html/freestat/rcomp/tmp/51mmv1293483851.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/html/freestat/rcomp/createtable file can be downloaded at http://www.wessa.net/cretab > load(file="/var/www/html/freestat/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/freestat/rcomp/tmp/6xw351293483851.tab") > > try(system("convert tmp/1fl671293483851.ps tmp/1fl671293483851.png",intern=TRUE)) character(0) > try(system("convert tmp/2qc5t1293483851.ps tmp/2qc5t1293483851.png",intern=TRUE)) character(0) > try(system("convert tmp/31mmv1293483851.ps tmp/31mmv1293483851.png",intern=TRUE)) character(0) > try(system("convert tmp/41mmv1293483851.ps tmp/41mmv1293483851.png",intern=TRUE)) character(0) > try(system("convert tmp/51mmv1293483851.ps tmp/51mmv1293483851.png",intern=TRUE)) character(0) > > > proc.time() user system elapsed 2.289 1.222 2.504