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(13328,12873,14000,13477,14237,13674,13529,14058,12975,14326,14008,16193,14483,14011,15057,14884,15414,14440,14900,15074,14442,15307,14938,17193,15528,14765,15838,15723,16150,15486,15986,15983,15692,16490,15686,18897,16316,15636,17163,16534,16518,16375,16290,16352,15943,16362,16393,19051,16747,16320,17910,16961,17480,17049,16879,17473,16998,17307,17418,20169,17871,17226,19062,17804,19100,18522,18060,18869,18127,18871,18890,21263,19547,18450,20254,19240,20216,19420,19415,20018,18652,19978,19509,21971) > par20 = '' > par19 = '' > par18 = '' > par17 = '' > par16 = '' > par15 = '' > par14 = '' > par13 = '' > par12 = '' > par11 = '' > par10 = '' > par9 = '' > par8 = '' > par7 = '' > par6 = '' > par5 = '' > par4 = '' > par3 = '' > par2 = '' > par1 = '12' > ylab = '' > xlab = '' > main = '' > #'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=52 ERROR: ABNORMAL_TERMINATION_IN_LNSRCH > m$coef level slope seas epsilon 15123.92748 51.98626 56600.95553 0.00000 > m$fitted level slope sea Jan 1 13328.00 0.000000 0.00000 Feb 1 13200.61 3.821522 -327.60944 Mar 1 13460.25 21.998810 539.74729 Apr 1 13535.99 26.277851 -58.98694 May 1 13789.38 42.162186 447.61931 Jun 1 13834.08 42.304814 -160.07855 Jul 1 13756.19 36.769578 -227.19417 Aug 1 13824.11 38.041365 233.89191 Sep 1 13582.50 27.030866 -607.49844 Oct 1 13729.75 31.841515 596.25292 Nov 1 13870.00 36.341621 138.00499 Dec 1 14657.75 68.754524 1535.24897 Jan 2 14940.97 73.376454 -457.96802 Feb 2 14973.30 72.310111 -962.30469 Mar 2 14910.20 66.739968 146.79647 Apr 2 14966.77 66.208983 -82.76597 May 2 14983.38 63.396575 430.61540 Jun 2 14841.90 51.759534 -401.90126 Jul 2 14877.57 50.877898 22.42794 Aug 2 14855.28 47.043254 218.71737 Sep 2 14996.14 51.770912 -554.13844 Oct 2 15086.87 53.672933 220.12612 Nov 2 15227.87 57.816811 -289.86755 Dec 2 15465.19 66.118613 1727.81212 Jan 3 15699.36 73.784964 -171.35685 Feb 3 15768.90 73.586433 -1003.89690 Mar 3 15770.23 70.022797 67.77490 Apr 3 15767.79 66.258555 -44.79401 May 3 15725.55 60.428441 424.45022 Jun 3 15780.16 60.111469 -294.15513 Jul 3 15877.23 62.117576 108.77136 Aug 3 15928.55 61.538676 54.45223 Sep 3 16100.19 67.347180 -408.19146 Oct 3 16281.84 73.282384 208.16444 Nov 3 16325.10 71.743292 -639.10099 Dec 3 16660.86 85.183051 2236.13798 Jan 4 16722.82 84.000666 -406.82047 Feb 4 16713.58 79.216302 -1077.57660 Mar 4 16831.78 81.242313 331.22473 Apr 4 16796.04 75.080656 -262.03581 May 4 16595.76 60.429461 -77.76321 Jun 4 16599.69 57.407487 -224.69347 Jul 4 16513.57 49.727294 -223.57406 Aug 4 16484.46 45.521619 -132.45883 Sep 4 16491.13 43.459935 -548.13240 Oct 4 16429.67 37.921934 -67.66548 Nov 4 16668.18 48.461666 -275.18054 Dec 4 16787.46 52.173488 2263.53739 Jan 5 16927.77 56.792678 -180.77229 Feb 5 17113.10 63.544176 -793.10417 Mar 5 17270.62 68.497715 639.37927 Apr 5 17247.87 63.669221 -286.86541 May 5 17325.07 64.387495 154.93197 Jun 5 17312.03 60.270933 -263.02777 Jul 5 17242.32 53.358151 -363.32330 Aug 5 17351.46 56.322600 121.53658 Sep 5 17478.67 60.084064 -480.67195 Oct 5 17568.38 61.653525 -261.38424 Nov 5 17693.24 64.996428 -275.23897 Dec 5 17860.16 70.382707 2308.84079 Jan 6 18029.55 75.614018 -158.54530 Feb 6 18110.44 75.893217 -884.44018 Mar 6 18234.47 78.440795 827.52991 Apr 6 18242.10 74.688995 -438.09918 May 6 18466.82 82.645578 633.17685 Jun 6 18650.02 87.981063 -128.02240 Jul 6 18689.53 85.408372 -629.52548 Aug 6 18776.60 85.496544 92.40438 Sep 6 18804.52 82.443164 -677.52057 Oct 6 18980.91 87.423264 -109.91304 Nov 6 19161.23 92.344869 -271.22632 Dec 6 19208.34 89.949116 2054.65999 Jan 7 19420.03 96.396648 126.97383 Feb 7 19505.23 95.803932 -1055.23322 Mar 7 19548.79 93.035766 705.20663 Apr 7 19679.91 95.053945 -439.90585 May 7 19740.47 93.225364 475.53014 Jun 7 19722.55 87.331864 -302.55414 Jul 7 19854.63 89.704438 -439.62596 Aug 7 19957.11 90.382277 60.88693 Sep 7 19856.40 80.250345 -1204.40444 Oct 7 19952.88 81.110564 25.11871 Nov 7 19964.53 77.428629 -455.52835 Dec 7 20006.57 75.553306 1964.42661 > m$resid Jan Feb Mar Apr May Jun 1 0.00000000 -1.30910910 1.31766719 0.31190190 1.58143678 0.01924012 2 1.70441854 -0.32901438 -1.01366319 -0.07218387 -0.35218043 -1.49386555 3 1.26383314 -0.03187679 -0.53598968 -0.52966684 -0.78914526 -0.04257228 4 -0.17201646 -0.69027057 0.28759722 -0.85920318 -2.01772362 -0.41444413 5 0.64952513 0.94730169 0.69187267 -0.67085706 0.09939941 -0.56879049 6 0.72843455 0.03885871 0.35411560 -0.52066382 1.10276968 0.73908206 7 0.89525039 -0.08229344 -0.38420072 0.27996970 -0.25356064 -0.81706395 Jul Aug Sep Oct Nov Dec 1 -0.93687478 0.24393278 -2.18563522 0.93563515 0.83989233 5.79778333 2 -0.11989363 -0.55070359 0.70732067 0.29317699 0.65522372 1.34613900 3 0.27249982 -0.08004306 0.81698551 0.84680444 -0.22206943 1.95350171 4 -1.05580443 -0.58126401 -0.28654054 -0.77335758 1.47741433 0.52164469 5 -0.95569707 0.41050292 0.52176269 0.21801503 0.46490634 0.74974825 6 -0.35642680 0.01222168 -0.42345570 0.69094569 0.68305951 -0.33259818 7 0.32894162 0.09399266 -1.40517313 0.11931565 -0.51074483 -0.26016137 > 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 > postscript(file="/var/www/html/freestat/rcomp/tmp/1q4d71291989353.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') > grid() > dev.off() null device 1 > mylagmax <- nx/2 > postscript(file="/var/www/html/freestat/rcomp/tmp/2q4d71291989353.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/3jvda1291989353.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/4unuc1291989353.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/5unuc1291989353.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/6qws31291989353.tab") > > try(system("convert tmp/1q4d71291989353.ps tmp/1q4d71291989353.png",intern=TRUE)) character(0) > try(system("convert tmp/2q4d71291989353.ps tmp/2q4d71291989353.png",intern=TRUE)) character(0) > try(system("convert tmp/3jvda1291989353.ps tmp/3jvda1291989353.png",intern=TRUE)) character(0) > try(system("convert tmp/4unuc1291989353.ps tmp/4unuc1291989353.png",intern=TRUE)) character(0) > try(system("convert tmp/5unuc1291989353.ps tmp/5unuc1291989353.png",intern=TRUE)) character(0) > > > proc.time() user system elapsed 2.607 1.161 2.769