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(1856,1834,2095,2164,2368,2072,2521,1823,1947,2226,1754,1786,2072,1846,2137,2466,2154,2289,2628,2074,2798,2194,2442,2565,2063,2069,2539,1898,2139,2408,2725,2201,2311,2548,2276,2351,2280,2057,2479,2379,2295,2456,2546,2844,2260,2981,2678,3440,2842,2450,2669,2570,2540,2318,2930,2947,2799,2695,2498,2260,2160,2058,2533,2150,2172,2155,3016,2333,2355,2825,2214,2360,2299,1746,2069,2267,1878,2266,2282,2085,2277,2251,1828,1954,1851,1570,1852,2187,1855,2218,2253,2028,2169,1997,2034,1791,1627,1631,2319,1707,1747,2397,2059,2251,2558,2406,2049,2074,1734) > 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') > m$coef level slope seas epsilon 5405.871899 5.795353 2663.357942 38347.350573 > m$fitted level slope sea Jan 1 1856.000 0.00000000 0.00000000 Feb 1 1849.007 -0.86627732 -2.24329947 Mar 1 1932.983 7.91315360 30.30376664 Apr 1 2016.690 14.20410734 25.10149358 May 1 2140.499 21.80000682 43.15130171 Jun 1 2130.159 19.86727950 -2.25688091 Jul 1 2261.333 25.86268097 61.19009551 Aug 1 2136.320 18.39740083 -39.42396557 Sep 1 2077.514 14.81954004 11.43276062 Oct 1 2126.387 16.31930689 36.43839553 Nov 1 2017.511 11.02274014 -29.77534765 Dec 1 1944.169 7.56725542 0.04496197 Jan 2 1990.990 4.20691014 -16.16710472 Feb 2 1954.188 2.45128101 -44.10907847 Mar 2 2009.062 4.91904677 46.49790702 Apr 2 2159.926 11.14489851 64.53378282 May 2 2163.905 10.86806298 2.46546560 Jun 2 2218.099 12.41887351 -5.75741568 Jul 2 2330.585 15.80985083 117.81150868 Aug 2 2281.139 13.68107396 -89.02729481 Sep 2 2439.350 18.27531375 95.74821970 Oct 2 2371.678 15.59101086 -20.87681978 Nov 2 2402.101 16.04687102 12.78585108 Dec 2 2464.697 17.36245896 14.70726128 Jan 3 2402.041 18.22624974 -177.71758741 Feb 3 2326.653 15.28968048 -99.19941374 Mar 3 2389.441 17.05521437 71.98232488 Apr 3 2234.486 10.88557746 -47.96339525 May 3 2203.888 9.48167372 6.39970240 Jun 3 2275.638 11.49051109 23.73655825 Jul 3 2381.092 14.42065724 178.56260123 Aug 3 2382.639 14.02818654 -158.87314809 Sep 3 2336.758 12.22834794 80.48214930 Oct 3 2409.982 14.04169093 29.67702805 Nov 3 2384.508 12.88707301 -38.21983316 Dec 3 2368.589 12.11906461 33.89231246 Jan 4 2385.423 12.17381233 -114.17503135 Feb 4 2327.943 10.16145841 -151.06528322 Mar 4 2336.668 10.11349803 144.72224827 Apr 4 2366.862 10.78269493 -21.64302913 May 4 2360.073 10.21482018 -35.06856918 Jun 4 2395.730 11.01316027 16.38057952 Jul 4 2395.611 10.67146609 169.70703092 Aug 4 2573.397 15.72476583 -20.37917103 Sep 4 2488.950 12.72624470 -54.20883096 Oct 4 2620.761 16.25895406 152.28765970 Nov 4 2665.736 17.09538428 -37.92501173 Dec 4 2893.817 22.83900824 176.52611023 Jan 5 2938.687 23.30392609 -135.97312321 Feb 5 2862.108 20.43660123 -239.89383149 Mar 5 2771.687 16.92546598 83.68285219 Apr 5 2720.361 14.74183909 -35.26558627 May 5 2685.029 13.16282549 -59.79673753 Jun 5 2581.092 9.53422685 -62.34218659 Jul 5 2645.216 11.20185694 190.79348977 Aug 5 2742.314 13.79885324 56.44272210 Sep 5 2804.958 15.26428230 -90.36182156 Oct 5 2753.511 13.27901722 56.83296881 Nov 5 2717.707 11.84119703 -134.82188404 Dec 5 2538.595 6.48111768 51.89478033 Jan 6 2444.554 3.91505880 -108.92760410 Feb 6 2381.705 1.97750226 -208.80409663 Mar 6 2387.042 2.08106807 140.27964660 Apr 6 2321.484 -0.02952577 -57.16385984 May 6 2280.315 -1.30568992 -38.37650835 Jun 6 2265.233 -1.72880042 -86.69289381 Jul 6 2432.573 3.41705771 293.59546172 Aug 6 2405.636 2.49979619 -20.51163959 Sep 6 2409.221 2.53239547 -56.08607667 Oct 6 2504.089 5.28732946 162.19872280 Nov 6 2448.348 3.48802115 -129.48395367 Dec 6 2395.471 1.86732456 61.38123302 Jan 7 2388.720 1.62873781 -74.82655966 Feb 7 2262.597 -2.12063781 -297.25685822 Mar 7 2146.951 -5.57380210 114.59150638 Apr 7 2180.801 -4.36260147 19.44581874 May 7 2105.384 -6.54164655 -106.63311377 Jun 7 2182.647 -3.98550420 -59.60530004 Jul 7 2127.306 -5.54307626 242.51425592 Aug 7 2116.220 -5.71034662 -21.73071551 Sep 7 2181.088 -3.58981717 -24.98150401 Oct 7 2149.468 -4.42757764 149.53561728 Nov 7 2082.939 -6.26835945 -148.65607163 Dec 7 2014.266 -8.09194155 46.51100156 Jan 8 1962.564 -9.34674069 -36.64275256 Feb 8 1913.950 -10.50865683 -276.67543417 Mar 8 1854.711 -11.98077700 80.08168503 Apr 8 1926.655 -9.42691386 118.05004035 May 8 1942.192 -8.66731000 -129.61498725 Jun 8 2029.566 -5.75339035 24.76377228 Jul 8 2028.008 -5.62653891 217.82861535 Aug 8 2036.938 -5.18788625 -33.81421009 Sep 8 2074.201 -3.91275451 22.23739445 Oct 8 2001.135 -5.98202639 114.01739906 Nov 8 2036.442 -4.75329286 -72.92968107 Dec 8 1946.622 -7.26467637 -10.42737755 Jan 9 1851.548 -9.84057110 -74.21572906 Feb 9 1852.157 -9.52971464 -239.02890505 Mar 9 1966.091 -5.81293480 142.95957916 Apr 9 1867.149 -8.63045763 -2.10266222 May 9 1867.950 -8.34496784 -136.97920958 Jun 9 2007.412 -3.87760145 137.85524743 Jul 9 1969.160 -4.91439061 148.47536514 Aug 9 2056.937 -2.12483963 35.84749371 Sep 9 2191.966 1.99342062 131.93326112 Oct 9 2233.480 3.17690303 105.10506727 Nov 9 2197.784 2.01692406 -82.52893370 Dec 9 2157.877 0.77139389 -12.42294043 Jan 10 2060.045 -2.15036893 -157.59600654 > m$resid Jan Feb Mar Apr May Jun 1 0.00000000 -0.08557131 0.92077376 0.84563664 1.26244009 -0.37984528 2 0.71004814 -0.45299505 0.58232971 1.70811593 -0.08677126 0.53509364 3 -1.16202972 -1.12764029 0.55906166 -2.07263574 -0.50987229 0.77463554 4 0.06339128 -0.86133732 -0.01732922 0.24487199 -0.21697365 0.31675483 5 0.28589944 -1.24575796 -1.35570839 -0.83806495 -0.61983686 -1.45797981 6 -1.28159821 -0.83476949 0.04138964 -0.83428851 -0.51010945 -0.17154139 7 -0.10884851 -1.59790282 -1.40517555 0.48774444 -0.88213129 1.04374403 8 -0.54801032 -0.49096836 -0.60477222 1.04044935 0.31020742 1.19633065 9 -1.10015618 0.13058117 1.53466733 -1.15618207 0.11727600 1.84135890 10 -1.23315431 Jul Aug Sep Oct Nov Dec 1 1.34105054 -1.84325537 -0.95255264 0.42307653 -1.56298931 -1.05681732 2 1.25006118 -0.82063689 1.82465241 -1.08751312 0.18793608 0.59178133 3 1.17701918 -0.16188814 -0.75498019 0.76956943 -0.49895240 -0.36472430 4 -0.13924737 2.09580208 -1.25797540 1.49648733 0.36096013 2.65674966 5 0.68198606 1.07500750 0.61182961 -0.83582676 -0.61484682 -2.39453214 6 2.11056718 -0.37939310 0.01357049 1.15451019 -0.76265453 -0.70492577 7 -0.64083016 -0.06921893 0.88158459 -0.34998675 -0.77487750 -0.77918521 8 0.05233914 0.18169785 0.52986283 -0.86266538 0.51471036 -1.06114207 9 -0.42875798 1.15660546 1.71102448 0.49269782 -0.48431281 -0.52269201 10 > 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/1t3271292507239.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/2t3271292507239.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/3muks1292507239.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/4fl1v1292507239.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/5fl1v1292507239.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/6tvzl1292507239.tab") > > try(system("convert tmp/1t3271292507239.ps tmp/1t3271292507239.png",intern=TRUE)) character(0) > try(system("convert tmp/2t3271292507239.ps tmp/2t3271292507239.png",intern=TRUE)) character(0) > try(system("convert tmp/3muks1292507239.ps tmp/3muks1292507239.png",intern=TRUE)) character(0) > try(system("convert tmp/4fl1v1292507239.ps tmp/4fl1v1292507239.png",intern=TRUE)) character(0) > try(system("convert tmp/5fl1v1292507239.ps tmp/5fl1v1292507239.png",intern=TRUE)) character(0) > > > proc.time() user system elapsed 2.831 1.191 3.009