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(6827,6178,7084,8162,8462,9644,10466,10748,9963,8194,6848,7027,7269,6775,7819,8371,9069,10248,11030,10882,10333,9109,7685,7602,8350,7829,8829,9948,10638,11253,11424,11391,10665,9396,7775,7933,8186,7444,8484,9864,10252,12282,11637,11577,12417,9637,8094,9280,8334,7899,9994,10078,10801,12950,12222,12246,13281,10366,8730,9614,8639,8772,10894,10455,11179,10588,10794,12770,13812,10857,9290,10925,9491,8919,11607,8852,12537,14759,13667,13731,15110,12185,10645,12161,10840,10436,13589,13402,13103,14933,14147,14057,16234,12389,11595,12772) > 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 285000.14 336028.51 45298.27 0.00 > m$fitted level slope sea Jan 1 6827.000 0.000000 0.000000 Feb 1 6235.730 -192.914149 -57.729993 Mar 1 6955.507 334.606297 128.493080 Apr 1 8120.634 826.267395 41.365974 May 1 8573.668 606.739518 -111.667934 Jun 1 9598.820 853.336396 45.180315 Jul 1 10503.976 883.903917 -37.975501 Aug 1 10852.220 567.889994 -104.220484 Sep 1 10145.320 -184.186557 -182.319519 Oct 1 8387.248 -1112.734872 -193.247987 Nov 1 6867.335 -1352.962029 -19.335074 Dec 1 6852.678 -563.383141 174.321510 Jan 2 6784.043 -277.873453 484.957484 Feb 2 6897.398 -65.426202 -122.397720 Mar 2 7688.125 412.435113 130.874577 Apr 2 8255.315 501.252960 115.685297 May 2 9181.880 744.146119 -112.879734 Jun 2 10220.845 912.561119 27.154948 Jul 2 11093.728 889.866830 -63.728286 Aug 2 11038.499 349.139748 -156.499385 Sep 2 10409.433 -210.544409 -76.432696 Oct 2 9296.553 -726.873493 -187.553311 Nov 2 7947.812 -1082.599537 -262.811695 Dec 2 7366.716 -796.580466 235.283524 Jan 3 7702.496 -149.710128 647.504382 Feb 3 8046.465 125.092140 -217.464628 Mar 3 8642.042 386.404228 186.958197 Apr 3 9792.204 819.166097 155.795547 May 3 10832.514 944.295975 -194.514332 Jun 3 11374.500 716.934735 -121.499585 Jul 3 11481.540 371.863159 -57.539738 Aug 3 11464.585 151.722907 -73.584586 Sep 3 10719.491 -356.136473 -54.491455 Oct 3 9536.194 -824.512272 -140.193614 Nov 3 8154.818 -1139.491117 -379.817608 Dec 3 7726.718 -737.896767 206.282264 Jan 4 7518.405 -437.948572 667.595433 Feb 4 7650.000 -119.630473 -205.999743 Mar 4 8335.618 327.981860 148.381593 Apr 4 9625.775 869.502748 238.225348 May 4 10436.662 836.482924 -184.661584 Jun 4 12171.623 1341.748418 110.376906 Jul 4 11982.534 480.288166 -345.534358 Aug 4 11622.526 7.005043 -45.526241 Sep 4 12164.071 308.168245 252.929347 Oct 4 10069.189 -1045.308897 -432.189410 Nov 4 8594.351 -1286.924866 -500.350646 Dec 4 8831.652 -430.054603 448.347771 Jan 5 7865.626 -732.023821 468.373500 Feb 5 8069.063 -209.173897 -170.062806 Mar 5 9722.778 827.291026 271.221714 Apr 5 10005.022 521.710753 72.978213 May 5 11109.972 849.321932 -308.971674 Jun 5 12477.300 1139.835849 472.699848 Jul 5 12657.349 601.458939 -435.348522 Aug 5 12553.829 205.678347 -307.829300 Sep 5 12515.588 68.693579 765.412172 Oct 5 10964.438 -840.614989 -598.437707 Nov 5 9565.776 -1153.541490 -835.776473 Dec 5 8888.036 -886.687745 725.963512 Jan 6 8306.691 -715.229886 332.309022 Feb 6 9064.776 107.659102 -292.776226 Mar 6 10315.118 743.712779 578.881927 Apr 6 10636.647 507.562290 -181.646983 May 6 11563.699 742.703208 -384.698931 Jun 6 10333.197 -362.048225 254.802668 Jul 6 10934.755 177.353239 -140.755063 Aug 6 12772.694 1107.527679 -2.694440 Sep 6 12882.624 548.552835 929.376009 Oct 6 11596.456 -479.027398 -739.455728 Nov 6 10281.999 -946.586658 -991.998968 Dec 6 10077.623 -530.981463 847.376589 Jan 7 9508.845 -552.157863 -17.845147 Feb 7 9392.139 -309.128010 -473.139013 Mar 7 10680.298 580.099475 926.701530 Apr 7 9446.297 -433.122714 -594.297323 May 7 11898.710 1181.785290 638.289570 Jun 7 14556.483 2007.170202 202.517245 Jul 7 14626.106 924.119556 -959.106041 Aug 7 13840.957 -31.742140 -109.957036 Sep 7 13765.214 -56.354417 1344.786082 Oct 7 12882.081 -518.647148 -697.081093 Nov 7 11838.568 -812.006833 -1193.568325 Dec 7 11209.617 -709.627729 951.382534 Jan 8 10860.045 -508.222056 -20.044900 Feb 8 11042.197 -123.175698 -606.197408 Mar 8 12040.549 501.160679 1548.451136 Apr 8 14376.824 1525.063594 -974.823535 May 8 13355.299 101.610023 -252.299373 Jun 8 14172.951 501.657290 760.048942 Jul 8 14779.805 560.398875 -632.804664 Aug 8 14355.881 10.644992 -298.880927 Sep 8 14626.337 155.774768 1607.663198 Oct 8 13220.235 -716.452102 -831.235448 Nov 8 12705.963 -603.560447 -1110.962710 Dec 8 11890.367 -722.030504 881.632782 > m$resid Jan Feb Mar Apr May Jun 1 0.00000000 -0.57477809 0.95389447 0.82330154 -0.37679723 0.42532080 2 0.54268818 0.34450788 0.84101860 0.15192315 0.41642292 0.29034966 3 1.14589954 0.46423163 0.45512130 0.74645916 0.21468110 -0.39171430 4 0.52222799 0.54388137 0.77564533 0.93652884 -0.05674597 0.87013707 5 -0.52277820 0.89711114 1.79205098 -0.52868662 0.56386405 0.50026511 6 0.29618499 1.41464505 1.09861881 -0.40844391 0.40514522 -1.90266618 7 -0.03654638 0.41822898 1.53513980 -1.75172504 2.78440660 1.42187860 8 0.34743018 0.66303216 1.07755125 1.76952409 -2.45536278 0.68933858 Jul Aug Sep Oct Nov Dec 1 0.05273226 -0.54515321 -1.29739965 -1.60182985 -0.41441363 1.36209520 2 -0.03914995 -0.93280073 -0.96550732 -0.89073497 -0.61365172 0.49450808 3 -0.59527589 -0.37975934 -0.87611229 -0.80801987 -0.54338684 0.69535545 4 -1.48593000 -0.81645744 0.51954186 -2.33487215 -0.41694328 1.48341156 5 -0.92848889 -0.68275735 -0.23631317 -1.56863612 -0.54015235 0.46166270 6 0.93010642 1.60458917 -0.96428233 -1.77275262 -0.80720883 0.71851544 7 -1.86731902 -1.64881360 -0.04245832 -0.79759797 -0.50650187 0.17689913 8 0.10127120 -0.94823687 0.25036269 -1.50499024 0.19491498 -0.20461626 > 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/14q0j1322579285.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/2wr4y1322579285.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/3663l1322579285.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/4dqwg1322579285.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/5nzr91322579285.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/60rkt1322579285.tab") > > try(system("convert tmp/14q0j1322579285.ps tmp/14q0j1322579285.png",intern=TRUE)) character(0) > try(system("convert tmp/2wr4y1322579285.ps tmp/2wr4y1322579285.png",intern=TRUE)) character(0) > try(system("convert tmp/3663l1322579285.ps tmp/3663l1322579285.png",intern=TRUE)) character(0) > try(system("convert tmp/4dqwg1322579285.ps tmp/4dqwg1322579285.png",intern=TRUE)) character(0) > try(system("convert tmp/5nzr91322579285.ps tmp/5nzr91322579285.png",intern=TRUE)) character(0) > > > proc.time() user system elapsed 2.720 0.372 3.064