R version 2.9.0 (2009-04-17) Copyright (C) 2009 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. 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(2540.9,2370.3,1807.5,1834.8,786.8,1561.4,1347.2,1549.8,1553.8,1822.5,3078.7,1589.1,1791.5,2558.1,2111.8,2083.1,2052.1,2243.5,2622,1952.6,808.9,1709.8,1582.1,865.6,1116.1,1119.4,2350,1975.6,2536.5,2785,2819.7,1829.5,758.3,2921.6,2482,1892.7,1855.1,2151.3,1642.2,1640.5,1366.1,1532.8,824.4,-518.7,-978.5,1162.5,1243.4,1199.5,883.1,1437.2,534.5,-1901.9,-2521.1,-1721.1,-3094.5,-3694.8,-2492.1,-464.6,-626.1,-1711.4) > par1 = '4' > #'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 412389.1366 279.1364 6777.4138 8891.4387 > m$fitted level slope sea 1 Q1 2540.9000 0.000000 0.000000 1 Q2 2403.2492 -33.026299 -32.810602 1 Q3 1899.0857 -60.114120 -82.234855 1 Q4 1882.6497 -58.316143 -48.731575 2 Q1 951.7291 29.483664 -147.250015 2 Q2 1434.8088 69.311405 120.441377 2 Q3 1332.3802 62.951185 18.101958 2 Q4 1432.8440 63.670298 116.223871 3 Q1 1761.9922 54.713415 -213.408600 3 Q2 1717.4225 49.370862 106.679397 3 Q3 2902.7012 86.837493 154.835762 3 Q4 1679.0502 67.194571 -64.958399 4 Q1 1949.6371 64.765058 -162.043353 4 Q2 2506.8468 84.263319 42.816323 4 Q3 1944.2956 64.329629 179.386024 4 Q4 2134.6302 66.198441 -53.897242 5 Q1 2263.7955 66.174772 -212.885108 5 Q2 2159.9025 60.653321 86.590907 5 Q3 2424.0447 66.583909 194.243876 5 Q4 2076.8602 59.963100 -116.564295 6 Q1 1144.7552 53.102423 -317.330582 6 Q2 1574.9931 63.903307 128.077680 6 Q3 1380.1496 56.703266 206.653722 6 Q4 941.3983 48.019846 -66.645733 7 Q1 1405.9602 52.970260 -297.574232 7 Q2 1029.3572 41.565779 97.763903 7 Q3 1951.6890 65.204468 382.324558 7 Q4 2087.1208 66.539791 -112.811878 8 Q1 2710.5445 75.120037 -184.297455 8 Q2 2794.4392 75.343144 -9.597546 8 Q3 2504.9415 65.807249 321.376837 8 Q4 2077.3582 55.781390 -238.819699 9 Q1 1089.7487 37.041526 -312.318146 9 Q2 2647.7745 74.852122 246.312313 9 Q3 2189.9407 61.180623 301.720827 9 Q4 2044.7606 56.762762 -148.290043 10 Q1 2310.7420 60.912676 -459.466680 10 Q2 1912.3441 49.608298 247.275398 10 Q3 1401.2984 35.385963 251.070393 10 Q4 1756.8269 42.521712 -122.165240 11 Q1 1772.7491 41.957143 -406.163857 11 Q2 1299.4052 29.309040 242.735406 11 Q3 702.5768 13.533094 133.180144 11 Q4 -259.9180 -8.906465 -241.008013 12 Q1 -599.2611 -16.262618 -373.220481 12 Q2 645.3026 14.725377 494.333733 12 Q3 958.4801 22.219967 279.506350 12 Q4 1356.0231 31.056599 -163.350666 13 Q1 1408.6091 31.552656 -525.900810 13 Q2 1041.3740 21.718500 403.059030 13 Q3 397.4228 5.019682 149.151787 13 Q4 -1461.0418 -39.664063 -406.986281 14 Q1 -1991.1192 -51.252102 -521.066134 14 Q2 -2197.5083 -55.093932 479.222395 14 Q3 -3305.7671 -81.525070 230.370331 14 Q4 -3359.7480 -80.855219 -335.552322 15 Q1 -2187.7493 -50.691881 -327.110698 15 Q2 -1212.0409 -25.162269 728.822442 15 Q3 -888.6204 -16.402476 256.197649 15 Q4 -1176.8404 -23.085444 -529.624050 > m$resid Qtr1 Qtr2 Qtr3 Qtr4 1 0.00000000 -0.05157600 -0.74644363 0.06919875 2 -1.54486562 0.59627521 -0.26844890 0.05861896 3 0.43362658 -0.14193812 1.76428977 -2.03750833 4 0.32321994 0.72699377 -1.00026155 0.19533926 5 0.09870193 -0.25513757 0.31374764 -0.64008558 6 -1.54300529 0.57083709 -0.39816349 -0.76511241 7 0.64479791 -0.65356960 1.35362416 0.10828835 8 0.85944570 0.01339105 -0.56026172 -0.75983968 9 -1.60708060 2.32546493 -0.81758643 -0.31747500 10 0.32181456 -0.70303018 -0.86021902 0.49211170 11 -0.04087637 -0.78925605 -0.96046653 -1.49932162 12 -0.50745352 1.93189343 0.45773530 0.57625174 13 0.03304684 -0.61117806 -1.02080702 -2.85994318 14 -0.75249705 -0.23779050 -1.61487216 0.04225950 15 1.92188281 1.57333725 0.53445710 -0.41693551 > 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/rcomp/tmp/14apy1293458085.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/rcomp/tmp/24apy1293458085.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/rcomp/tmp/3x17j1293458085.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/rcomp/tmp/4x17j1293458085.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/rcomp/tmp/5psom1293458085.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/rcomp/createtable file can be downloaded at http://www.wessa.net/cretab > load(file="/var/www/html/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/rcomp/tmp/63k4v1293458085.tab") > > try(system("convert tmp/14apy1293458085.ps tmp/14apy1293458085.png",intern=TRUE)) character(0) > try(system("convert tmp/24apy1293458085.ps tmp/24apy1293458085.png",intern=TRUE)) character(0) > try(system("convert tmp/3x17j1293458085.ps tmp/3x17j1293458085.png",intern=TRUE)) character(0) > try(system("convert tmp/4x17j1293458085.ps tmp/4x17j1293458085.png",intern=TRUE)) character(0) > try(system("convert tmp/5psom1293458085.ps tmp/5psom1293458085.png",intern=TRUE)) character(0) > > > proc.time() user system elapsed 1.207 0.823 3.094