R version 2.12.1 (2010-12-16) Copyright (C) 2010 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(547612,563280,581302,572273,518654,520579,530577,540324,547970,555654,551174,548604,563668,586111,604378,600991,544686,537034,551531,563250,574761,580112,575093,557560,564478,580523,596594,586570,536214,523597,536535,536322,532638,528222,516141,501866,506174,517945,533590,528379,477580,469357,490243,492622,507561,516922,514258,509846,527070,541657,564591,555362,498662,511038,525919,531673,548854,560576,557274,565742,587625,619916,625809,619567,572942,572775,574205,579799,590072,593408,597141,595404,612117,628232,628884,620735,569028,567456,573100,584428,589379,590865,595454,594167,611324,612613,610763,593530,542722,536662,543599,555332,560854,562325,554788,547344,565464,577992,579714,569323,506971,500857,509127,509933,517009,519164,512238,509239,518585,522975,525192,516847,455626,454724,461251,470439,474605,476049,471067,470984,502831,512927,509673,484015,431328,436087,442867,447988,460070,467037,460170,464196,485025) > 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 36831935.2 708401.7 35747328.4 0.0 > m$fitted level slope sea Jan 1 547612.0 0.0000 0.0000 Feb 1 557514.0 429.0194 5766.0113 Mar 1 572723.0 1497.6594 8579.0088 Apr 1 574844.9 1542.7214 -2571.8597 May 1 547297.5 -534.6611 -28643.4871 Jun 1 527087.5 -2033.3942 -6508.4569 Jul 1 523431.8 -2169.6065 7145.1820 Aug 1 530435.5 -1326.5588 9888.5226 Sep 1 540451.4 -208.1605 7518.6207 Oct 1 550240.9 829.8126 5413.0828 Nov 1 552864.6 1023.1569 -1690.5509 Dec 1 551204.2 725.8558 -2600.2391 Jan 2 561784.8 1716.5764 1883.2284 Feb 2 577640.1 3289.1679 8470.8962 Mar 2 589407.2 4283.4100 14970.8413 Apr 2 591464.8 4020.5371 9526.1553 May 2 579106.2 2090.8558 -34420.2475 Jun 2 559972.9 -404.4753 -22938.9192 Jul 2 551442.0 -1359.5819 88.9982 Aug 2 553657.2 -938.9071 9592.8060 Sep 2 562830.0 253.4100 11930.9954 Oct 2 570699.0 1152.6662 9412.9745 Nov 2 575907.7 1631.3323 -814.7342 Dec 2 573521.9 1158.5965 -15961.9283 Jan 3 572618.6 916.2952 -8140.5887 Feb 3 574372.4 1015.3538 6150.5764 Mar 3 576178.5 1109.1511 20415.4995 Apr 3 570753.2 335.4264 15816.7756 May 3 563481.0 -564.4525 -27266.9973 Jun 3 551837.1 -1875.9013 -28240.1194 Jul 3 543404.0 -2652.5899 -6869.0466 Aug 3 536077.6 -3206.1813 244.4100 Sep 3 528294.9 -3747.9548 4343.1472 Oct 3 521621.3 -4093.9880 6600.7459 Nov 3 515687.8 -4311.2963 453.1593 Dec 3 514042.2 -3996.6246 -12176.2160 Jan 4 512895.5 -3659.9820 -6721.4862 Feb 4 510530.1 -3506.8232 7414.9094 Mar 4 507839.7 -3410.2385 25750.2824 Apr 4 505681.3 -3262.3868 22697.6729 May 4 500392.7 -3501.4702 -22812.6630 Jun 4 494775.1 -3751.2582 -25418.1131 Jul 4 492043.2 -3630.8233 -1800.2204 Aug 4 489013.6 -3559.7650 3608.3851 Sep 4 493526.7 -2605.9557 14034.3260 Oct 4 501706.5 -1332.6277 15215.5452 Nov 4 509704.1 -232.0083 4553.8660 Dec 4 517983.9 771.8090 -8137.8887 Jan 5 527879.9 1848.4635 -809.9146 Feb 5 533474.1 2290.7692 8182.8836 Mar 5 537721.4 2521.7590 26869.5852 Apr 5 535772.6 1994.4674 19589.4458 May 5 528265.0 874.5400 -29602.9562 Jun 5 531736.5 1180.6636 -20698.4611 Jul 5 531820.6 1051.3249 -5901.5532 Aug 5 533565.9 1133.2136 -1892.8977 Sep 5 538765.0 1612.8638 10089.0435 Oct 5 546991.7 2392.6531 13584.3345 Nov 5 555270.1 3086.2264 2003.9035 Dec 5 569749.8 4428.6980 -4007.8309 Jan 6 583927.1 5577.7958 3697.9082 Feb 6 601715.6 7017.6881 18200.4299 Mar 6 603390.6 6387.7435 22418.4420 Apr 6 601566.5 5420.0595 18000.5278 May 6 604005.5 5068.9345 -31063.4529 Jun 6 600646.5 4076.2400 -27871.5048 Jul 6 591842.6 2558.5697 -17637.5766 Aug 6 588025.6 1807.1373 -8226.6323 Sep 6 586696.7 1437.5478 3375.3030 Oct 6 586612.0 1258.2188 6795.9958 Nov 6 595686.5 2178.6907 1454.5440 Dec 6 603523.5 2845.0120 -8119.5205 Jan 7 610729.0 3358.6245 1387.9814 Feb 7 610450.2 2930.0670 17781.8497 Mar 7 606701.6 2143.2649 22182.3687 Apr 7 602863.7 1438.8578 17871.3107 May 7 597701.5 661.6871 -28673.5153 Jun 7 591834.9 -106.9044 -24378.9497 Jul 7 588366.7 -502.7308 -15266.7091 Aug 7 589166.4 -349.3307 -4738.4043 Sep 7 588314.8 -408.4904 1064.2384 Oct 7 588761.6 -307.7785 2103.3988 Nov 7 593280.1 260.3725 2173.8579 Dec 7 599286.5 936.7703 -5119.5488 Jan 8 604126.2 1396.2733 7197.7847 Feb 8 598141.3 527.1271 14471.6838 Mar 8 590047.5 -487.9737 20715.4934 Apr 8 578890.7 -1743.9585 14639.2829 May 8 570764.6 -2495.1602 -28042.5675 Jun 8 563080.6 -3105.8908 -26418.5509 Jul 8 559131.8 -3205.1103 -15532.7848 Aug 8 558169.4 -2941.0802 -2837.3616 Sep 8 558906.6 -2508.0420 1947.3847 Oct 8 560917.2 -1976.1529 1407.7593 Nov 8 558181.4 -2065.5649 -3393.3773 Dec 8 554204.4 -2290.5133 -6860.4431 Jan 9 551721.8 -2313.1232 13742.1673 Feb 9 552921.9 -1899.5845 25070.1461 Mar 9 551917.5 -1794.2119 27796.5449 Apr 9 550260.2 -1778.0970 19062.8379 May 9 540351.0 -2734.9819 -33380.0078 Jun 9 531372.2 -3469.7632 -30515.1644 Jul 9 526534.4 -3630.7648 -17407.3806 Aug 9 519118.8 -4076.2263 -9185.8370 Sep 9 515620.3 -4008.2393 1388.7453 Oct 9 514492.4 -3669.2626 4671.6492 Nov 9 513313.4 -3376.2113 -1075.4123 Dec 9 513759.0 -2926.4814 -4520.0307 Jan 10 509385.3 -3096.7958 9199.7114 Feb 10 501849.1 -3619.2522 21125.8593 Mar 10 496201.7 -3857.9413 28990.2818 Apr 10 491990.3 -3899.5383 24856.7028 May 10 486637.7 -4070.5165 -31011.7312 Jun 10 482990.4 -4020.7197 -28266.4098 Jul 10 478068.2 -4126.8044 -16817.1747 Aug 10 477215.0 -3741.5681 -6775.9836 Sep 10 474968.4 -3565.6353 -363.4404 Oct 10 472575.3 -3427.6676 3473.7242 Nov 10 471519.8 -3148.5367 -452.8383 Dec 10 471375.5 -2795.0432 -391.5487 Jan 11 480116.9 -1437.5718 22714.0874 Feb 11 486229.1 -549.1647 26697.8503 Mar 11 483849.2 -764.5932 25823.7595 Apr 11 470300.4 -2268.8840 13714.5731 May 11 463348.3 -2819.9239 -32020.3220 Jun 11 461951.0 -2652.5401 -25863.9866 Jul 11 460769.9 -2479.3994 -17902.9258 Aug 11 457309.7 -2594.8066 -9321.7381 Sep 11 458123.6 -2193.7117 1946.3610 Oct 11 461590.0 -1527.7298 5447.0423 Nov 11 463335.5 -1142.5998 -3165.4737 Dec 11 467477.9 -520.7680 -3281.9108 Jan 12 466723.0 -548.3170 18301.9963 > m$resid Jan Feb Mar Apr May Jun 1 0.00000000 1.23416891 1.80689080 0.09276974 -4.49756231 -3.00861046 2 1.49576596 2.08513715 1.15148073 -0.30652400 -2.32284093 -3.02667588 3 -0.29532645 0.11959423 0.11035849 -0.90773832 -1.06847717 -1.56674750 4 0.40352855 0.18294606 0.11433619 0.17462333 -0.28385950 -0.29793245 5 1.28524601 0.52689182 0.27395706 -0.62444424 -1.32972635 0.36460052 6 1.36909041 1.71332001 -0.74772242 -1.14751513 -0.41695529 -1.18114937 7 0.61127255 -0.50963966 -0.93430285 -0.83594989 -0.92301000 -0.91394772 8 0.54651241 -1.03323564 -1.20568960 -1.49121166 -0.89227461 -0.72597294 9 -0.02688036 0.49150972 0.12517382 0.01913823 -1.13668504 -0.87325040 10 -0.20243183 -0.62088172 -0.28356480 -0.04940954 -0.20311777 0.05917386 11 1.61321393 1.05568570 -0.25594213 -1.78700017 -0.65464923 0.19888963 12 -0.03273621 Jul Aug Sep Oct Nov Dec 1 -0.24459284 1.36772408 1.67636960 1.46744392 0.26185700 -0.39002046 2 -1.15402455 0.50607741 1.42995599 1.07607965 0.57225715 -0.56774291 3 -0.92648360 -0.65822312 -0.64281206 -0.41027529 -0.25800949 0.37536501 4 0.14363888 0.08451030 1.13132589 1.50898338 1.30644736 1.19566297 5 -0.15412546 0.09739799 0.56920711 0.92465811 0.82347711 1.59760628 6 -1.80680834 -0.89356735 -0.43879735 -0.21278009 1.09310238 0.79246669 7 -0.47089633 0.18237122 -0.07025911 0.11955536 0.67480789 0.80412801 8 -0.11797930 0.31383164 0.51437986 0.63161099 -0.10620850 -0.26736061 9 -0.19138331 -0.52940799 0.08076628 0.40261324 0.34812809 0.53444289 10 -0.12607904 0.45778834 0.20901476 0.16388951 0.33160833 0.42004198 11 0.20574911 -0.13713310 0.47653135 0.79117126 0.45755127 0.73885989 12 > 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/rcomp/tmp/1td6s1322679264.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/rcomp/tmp/2avm91322679264.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/rcomp/tmp/3pniw1322679264.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/rcomp/tmp/4ocls1322679264.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/rcomp/tmp/5wtkd1322679264.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/rcomp/createtable file can be downloaded at http://www.wessa.net/cretab > load(file="/var/www/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/rcomp/tmp/6ovry1322679264.tab") > > try(system("convert tmp/1td6s1322679264.ps tmp/1td6s1322679264.png",intern=TRUE)) character(0) > try(system("convert tmp/2avm91322679264.ps tmp/2avm91322679264.png",intern=TRUE)) character(0) > try(system("convert tmp/3pniw1322679264.ps tmp/3pniw1322679264.png",intern=TRUE)) character(0) > try(system("convert tmp/4ocls1322679264.ps tmp/4ocls1322679264.png",intern=TRUE)) character(0) > try(system("convert tmp/5wtkd1322679264.ps tmp/5wtkd1322679264.png",intern=TRUE)) character(0) > > > proc.time() user system elapsed 3.272 0.356 3.616