R version 2.15.2 (2012-10-26) -- "Trick or Treat" Copyright (C) 2012 The R Foundation for Statistical Computing ISBN 3-900051-07-0 Platform: i686-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(9700,9081,9084,9743,8587,9731,9563,9998,9437,10038,9918,9252,9737,9035,9133,9487,8700,9627,8947,9283,8829,9947,9628,9318,9605,8640,9214,9567,8547,9185,9470,9123,9278,10170,9434,9655,9429,8739,9552,9687,9019,9672,9206,9069,9788,10312,10105,9863,9656,9295,9946,9701,9049,10190,9706,9765,9893,9994,10433,10073,10112,9266,9820,10097,9115,10411,9678,10408,10153,10368,10581,10597,10680,9738,9556) > par1 = '12' > 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 17714.3239 0.0000 34613.9280 808.0046 > m$fitted level slope sea Jan 1 9700.000 0.00000000 0.0000000 Feb 1 9420.726 -4.57650591 -334.7706850 Mar 1 9148.229 -22.22438754 -61.4635790 Apr 1 9343.354 -8.81461753 396.7847885 May 1 9095.714 -19.39029604 -504.9206510 Jun 1 9271.821 -13.67981790 455.8577752 Jul 1 9434.357 -10.03481868 125.5925470 Aug 1 9698.073 -5.33029079 295.1675501 Sep 1 9674.895 -5.62166050 -237.5848310 Oct 1 9807.965 -3.34401391 227.6276516 Nov 1 9889.115 -1.93839486 27.4196920 Dec 1 9667.907 -5.59764315 -412.1074455 Jan 2 9552.084 -3.60454864 186.8839290 Feb 2 9409.587 -2.67974580 -372.1245401 Mar 2 9311.853 -4.16578472 -177.3463537 Apr 2 9154.695 -8.33829474 334.5980597 May 2 9193.577 -6.97830915 -494.3022491 Jun 2 9253.457 -5.30776283 372.4747554 Jul 2 9171.221 -6.85029545 -222.9605002 Aug 2 9097.546 -7.92883151 186.5624231 Sep 2 9105.004 -7.72121733 -276.2599725 Oct 2 9330.421 -5.08184193 612.6854073 Nov 2 9424.025 -4.23497530 202.3253543 Dec 2 9518.242 -3.76506233 -201.8928656 Jan 3 9459.674 -3.82937373 146.2497889 Feb 3 9244.035 -4.56264528 -600.4958281 Mar 3 9214.613 -4.81131269 -0.2101678 Apr 3 9224.507 -4.57503226 342.2610226 May 3 9155.485 -5.80637816 -607.4694299 Jun 3 8982.933 -8.99054793 204.7299860 Jul 3 9195.957 -5.17696138 270.4387455 Aug 3 9183.677 -5.28013568 -60.5604467 Sep 3 9382.207 -2.87427936 -107.5756489 Oct 3 9524.112 -1.54142129 643.4880955 Nov 3 9459.841 -1.96807599 -24.7993399 Dec 3 9554.227 -1.49285886 99.1717364 Jan 4 9426.440 -2.04248226 4.6568155 Feb 4 9354.034 -2.44846525 -613.8746734 Mar 4 9400.483 -2.02148594 150.7213060 Apr 4 9364.459 -2.42573642 323.0874837 May 4 9439.233 -1.34211677 -421.4686772 Jun 4 9526.304 -0.04386046 144.2742107 Jul 4 9354.355 -2.44659442 -145.5658182 Aug 4 9291.780 -3.19071120 -221.7961331 Sep 4 9527.285 -0.71668869 256.7889156 Oct 4 9624.440 0.09932692 685.9447506 Nov 4 9836.145 1.50094826 265.3580333 Dec 4 9824.456 1.42770128 38.7616635 Jan 5 9747.286 1.00155298 -89.9896545 Feb 5 9801.811 1.33980742 -507.6894910 Mar 5 9816.196 1.44365870 129.5912922 Apr 5 9676.986 0.07026656 26.2914908 May 5 9569.022 -1.13320563 -518.2771269 Jun 5 9686.250 0.25460997 501.8344412 Jul 5 9795.079 1.49863032 -90.8441274 Aug 5 9964.849 3.26699000 -202.5986916 Sep 5 9928.379 2.90252605 -34.7268605 Oct 5 9736.872 1.39225478 260.3269608 Nov 5 9841.188 2.07126977 590.1166885 Dec 5 9924.941 2.55338459 146.7131006 Jan 6 10052.813 3.28399086 57.1252793 Feb 6 9964.622 2.70104439 -697.1218961 Mar 6 9805.671 1.51386933 16.9685498 Apr 6 9863.903 1.99201048 232.1745911 May 6 9841.325 1.76344638 -725.9261269 Jun 6 9901.651 2.33314215 508.3969426 Jul 6 9895.089 2.24753785 -216.9441359 Aug 6 10131.531 4.36822758 272.6403054 Sep 6 10180.012 4.72970005 -27.7355548 Oct 6 10200.941 4.84730993 166.7926682 Nov 6 10156.271 4.52693611 425.5433138 Dec 6 10266.994 5.16344572 328.2593794 Jan 7 10388.595 5.85285141 289.4919777 Feb 7 10404.405 5.91494618 -666.5684505 Mar 7 10104.973 3.83205984 -543.9786813 > m$resid Jan Feb Mar Apr May Jun 1 0.00000000 -1.94832687 -1.35926283 1.35805753 -1.70074937 1.45002005 2 -0.86793921 -1.07955414 -0.67260263 -1.05056490 0.33340973 0.48735258 3 -0.41306202 -1.58338540 -0.18131278 0.10544201 -0.46351598 -1.21607829 4 -0.94398281 -0.52262084 0.35935855 -0.24786889 0.56251191 0.64812329 5 -0.58587764 0.39727801 0.09628608 -1.03368271 -0.79340520 0.87172896 6 0.93323451 -0.67934292 -1.19660502 0.41881886 -0.18132785 0.43282211 7 0.86689576 0.07400514 -2.26496731 Jul Aug Sep Oct Nov Dec 1 1.31901460 2.05145162 -0.13362994 1.03724070 0.63140501 -1.63795976 2 -0.57079721 -0.49948250 0.11528118 1.74695556 0.73898707 0.73815315 3 1.64065883 -0.05291751 1.52370115 1.08336130 -0.46942512 0.72113612 4 -1.27009465 -0.44696474 1.78035553 0.73100445 1.58103251 -0.09850993 5 0.80322135 1.24992982 -0.29596626 -1.44999649 0.76798107 0.60921981 6 -0.06591191 1.73986641 0.32837140 0.12072339 -0.36917555 0.79149662 7 > 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/fisher/rcomp/tmp/1ynkm1355494886.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/fisher/rcomp/tmp/23l8h1355494886.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/fisher/rcomp/tmp/332fb1355494886.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/fisher/rcomp/tmp/4bcob1355494886.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/fisher/rcomp/tmp/5n9w51355494886.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/fisher/rcomp/createtable file can be downloaded at http://www.wessa.net/cretab > load(file="/var/fisher/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/fisher/rcomp/tmp/6ld0g1355494886.tab") > > try(system("convert tmp/1ynkm1355494886.ps tmp/1ynkm1355494886.png",intern=TRUE)) character(0) > try(system("convert tmp/23l8h1355494886.ps tmp/23l8h1355494886.png",intern=TRUE)) character(0) > try(system("convert tmp/332fb1355494886.ps tmp/332fb1355494886.png",intern=TRUE)) character(0) > try(system("convert tmp/4bcob1355494886.ps tmp/4bcob1355494886.png",intern=TRUE)) character(0) > try(system("convert tmp/5n9w51355494886.ps tmp/5n9w51355494886.png",intern=TRUE)) character(0) > > > proc.time() user system elapsed 3.312 0.767 4.105