R version 2.12.0 (2010-10-15) 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 36831934.7 708401.7 35747328.2 0.0 > m$fitted level slope sea Jan 1 547612.0 0.0000 0.00000 Feb 1 557514.0 429.0194 5766.01137 Mar 1 572723.0 1497.6594 8579.00877 Apr 1 574844.9 1542.7214 -2571.85971 May 1 547297.5 -534.6611 -28643.48718 Jun 1 527087.5 -2033.3943 -6508.45695 Jul 1 523431.8 -2169.6066 7145.18204 Aug 1 530435.5 -1326.5588 9888.52268 Sep 1 540451.4 -208.1605 7518.62071 Oct 1 550240.9 829.8126 5413.08286 Nov 1 552864.6 1023.1570 -1690.55094 Dec 1 551204.2 725.8558 -2600.23907 Jan 2 561784.8 1716.5764 1883.22839 Feb 2 577640.1 3289.1679 8470.89624 Mar 2 589407.2 4283.4100 14970.84129 Apr 2 591464.8 4020.5371 9526.15526 May 2 579106.2 2090.8558 -34420.24755 Jun 2 559972.9 -404.4753 -22938.91924 Jul 2 551442.0 -1359.5819 88.99816 Aug 2 553657.2 -938.9072 9592.80601 Sep 2 562830.0 253.4100 11930.99538 Oct 2 570699.0 1152.6662 9412.97448 Nov 2 575907.7 1631.3323 -814.73423 Dec 2 573521.9 1158.5966 -15961.92828 Jan 3 572618.6 916.2952 -8140.58869 Feb 3 574372.4 1015.3538 6150.57642 Mar 3 576178.5 1109.1511 20415.49956 Apr 3 570753.2 335.4264 15816.77560 May 3 563481.0 -564.4525 -27266.99733 Jun 3 551837.1 -1875.9013 -28240.11946 Jul 3 543404.0 -2652.5899 -6869.04664 Aug 3 536077.6 -3206.1813 244.41002 Sep 3 528294.9 -3747.9549 4343.14717 Oct 3 521621.3 -4093.9880 6600.74591 Nov 3 515687.8 -4311.2963 453.15936 Dec 3 514042.2 -3996.6247 -12176.21595 Jan 4 512895.5 -3659.9821 -6721.48617 Feb 4 510530.1 -3506.8232 7414.90947 Mar 4 507839.7 -3410.2385 25750.28241 Apr 4 505681.3 -3262.3868 22697.67298 May 4 500392.7 -3501.4702 -22812.66295 Jun 4 494775.1 -3751.2582 -25418.11312 Jul 4 492043.2 -3630.8233 -1800.22041 Aug 4 489013.6 -3559.7650 3608.38509 Sep 4 493526.7 -2605.9557 14034.32601 Oct 4 501706.5 -1332.6276 15215.54518 Nov 4 509704.1 -232.0082 4553.86596 Dec 4 517983.9 771.8091 -8137.88875 Jan 5 527879.9 1848.4635 -809.91462 Feb 5 533474.1 2290.7692 8182.88355 Mar 5 537721.4 2521.7591 26869.58521 Apr 5 535772.6 1994.4674 19589.44577 May 5 528265.0 874.5400 -29602.95627 Jun 5 531736.5 1180.6636 -20698.46115 Jul 5 531820.6 1051.3249 -5901.55323 Aug 5 533565.9 1133.2137 -1892.89777 Sep 5 538765.0 1612.8638 10089.04351 Oct 5 546991.7 2392.6532 13584.33452 Nov 5 555270.1 3086.2264 2003.90351 Dec 5 569749.8 4428.6981 -4007.83087 Jan 6 583927.1 5577.7958 3697.90817 Feb 6 601715.6 7017.6882 18200.42985 Mar 6 603390.6 6387.7435 22418.44194 Apr 6 601566.5 5420.0596 18000.52775 May 6 604005.5 5068.9346 -31063.45295 Jun 6 600646.5 4076.2400 -27871.50488 Jul 6 591842.6 2558.5696 -17637.57661 Aug 6 588025.6 1807.1373 -8226.63238 Sep 6 586696.7 1437.5478 3375.30301 Oct 6 586612.0 1258.2188 6795.99582 Nov 6 595686.5 2178.6907 1454.54401 Dec 6 603523.5 2845.0120 -8119.52046 Jan 7 610729.0 3358.6245 1387.98138 Feb 7 610450.2 2930.0670 17781.84977 Mar 7 606701.6 2143.2648 22182.36877 Apr 7 602863.7 1438.8577 17871.31068 May 7 597701.5 661.6871 -28673.51524 Jun 7 591834.9 -106.9044 -24378.94972 Jul 7 588366.7 -502.7308 -15266.70908 Aug 7 589166.4 -349.3308 -4738.40432 Sep 7 588314.8 -408.4904 1064.23838 Oct 7 588761.6 -307.7786 2103.39878 Nov 7 593280.1 260.3725 2173.85795 Dec 7 599286.5 936.7703 -5119.54877 Jan 8 604126.2 1396.2733 7197.78473 Feb 8 598141.3 527.1271 14471.68378 Mar 8 590047.5 -487.9737 20715.49345 Apr 8 578890.7 -1743.9586 14639.28295 May 8 570764.6 -2495.1603 -28042.56744 Jun 8 563080.6 -3105.8909 -26418.55089 Jul 8 559131.8 -3205.1103 -15532.78476 Aug 8 558169.4 -2941.0802 -2837.36165 Sep 8 558906.6 -2508.0420 1947.38473 Oct 8 560917.2 -1976.1529 1407.75925 Nov 8 558181.4 -2065.5650 -3393.37733 Dec 8 554204.4 -2290.5133 -6860.44305 Jan 9 551721.8 -2313.1232 13742.16736 Feb 9 552921.9 -1899.5845 25070.14614 Mar 9 551917.5 -1794.2119 27796.54490 Apr 9 550260.2 -1778.0971 19062.83789 May 9 540351.0 -2734.9819 -33380.00776 Jun 9 531372.2 -3469.7633 -30515.16445 Jul 9 526534.4 -3630.7648 -17407.38059 Aug 9 519118.8 -4076.2264 -9185.83701 Sep 9 515620.3 -4008.2393 1388.74526 Oct 9 514492.4 -3669.2626 4671.64919 Nov 9 513313.4 -3376.2113 -1075.41232 Dec 9 513759.0 -2926.4814 -4520.03068 Jan 10 509385.3 -3096.7959 9199.71145 Feb 10 501849.1 -3619.2522 21125.85933 Mar 10 496201.7 -3857.9413 28990.28179 Apr 10 491990.3 -3899.5384 24856.70287 May 10 486637.7 -4070.5165 -31011.73121 Jun 10 482990.4 -4020.7197 -28266.40978 Jul 10 478068.2 -4126.8044 -16817.17476 Aug 10 477215.0 -3741.5681 -6775.98359 Sep 10 474968.4 -3565.6353 -363.44046 Oct 10 472575.3 -3427.6676 3473.72414 Nov 10 471519.8 -3148.5367 -452.83827 Dec 10 471375.5 -2795.0432 -391.54868 Jan 11 480116.9 -1437.5718 22714.08741 Feb 11 486229.1 -549.1647 26697.85026 Mar 11 483849.2 -764.5932 25823.75947 Apr 11 470300.4 -2268.8840 13714.57310 May 11 463348.3 -2819.9239 -32020.32204 Jun 11 461951.0 -2652.5401 -25863.98659 Jul 11 460769.9 -2479.3994 -17902.92579 Aug 11 457309.7 -2594.8066 -9321.73809 Sep 11 458123.6 -2193.7117 1946.36103 Oct 11 461590.0 -1527.7298 5447.04231 Nov 11 463335.5 -1142.5998 -3165.47374 Dec 11 467477.9 -520.7679 -3281.91077 Jan 12 466723.0 -548.3170 18301.99627 > m$resid Jan Feb Mar Apr May Jun 1 0.00000000 1.23416892 1.80689081 0.09276975 -4.49756232 -3.00861048 2 1.49576596 2.08513715 1.15148072 -0.30652400 -2.32284094 -3.02667590 3 -0.29532645 0.11959422 0.11035849 -0.90773832 -1.06847717 -1.56674750 4 0.40352856 0.18294606 0.11433619 0.17462334 -0.28385950 -0.29793245 5 1.28524601 0.52689181 0.27395705 -0.62444425 -1.32972635 0.36460052 6 1.36909041 1.71332001 -0.74772243 -1.14751514 -0.41695529 -1.18114938 7 0.61127255 -0.50963967 -0.93430286 -0.83594989 -0.92301000 -0.91394772 8 0.54651241 -1.03323565 -1.20568960 -1.49121166 -0.89227461 -0.72597293 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.61321394 1.05568570 -0.25594214 -1.78700018 -0.65464924 0.19888964 12 -0.03273622 Jul Aug Sep Oct Nov Dec 1 -0.24459284 1.36772410 1.67636961 1.46744393 0.26185700 -0.39002046 2 -1.15402456 0.50607741 1.42995600 1.07607965 0.57225716 -0.56774292 3 -0.92648360 -0.65822312 -0.64281206 -0.41027529 -0.25800948 0.37536502 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.80680835 -0.89356735 -0.43879735 -0.21278009 1.09310238 0.79246669 7 -0.47089632 0.18237122 -0.07025911 0.11955536 0.67480789 0.80412801 8 -0.11797929 0.31383165 0.51437987 0.63161100 -0.10620850 -0.26736061 9 -0.19138331 -0.52940798 0.08076628 0.40261324 0.34812809 0.53444289 10 -0.12607904 0.45778835 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/1chz71322658019.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/2crjq1322658019.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/35z0d1322658019.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/4d4xj1322658019.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/5s6o51322658019.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/6b2dj1322658019.tab") > > try(system("convert tmp/1chz71322658019.ps tmp/1chz71322658019.png",intern=TRUE)) character(0) > try(system("convert tmp/2crjq1322658019.ps tmp/2crjq1322658019.png",intern=TRUE)) character(0) > try(system("convert tmp/35z0d1322658019.ps tmp/35z0d1322658019.png",intern=TRUE)) character(0) > try(system("convert tmp/4d4xj1322658019.ps tmp/4d4xj1322658019.png",intern=TRUE)) character(0) > try(system("convert tmp/5s6o51322658019.ps tmp/5s6o51322658019.png",intern=TRUE)) character(0) > > > proc.time() user system elapsed 2.650 0.170 2.807