R version 2.13.0 (2011-04-13) Copyright (C) 2011 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(183.046,175.714,175.768,171.029,170.465,170.102,156.389,124.291,99.360,86.675,85.056,128.236,164.257,162.401,152.779,156.005,153.387,153.190,148.840,144.211,145.953,145.542,150.271,147.489,143.824,134.754,131.736,126.304,125.511,125.495,130.133,126.257,110.323,98.417,105.749,120.665,124.075,127.245,146.731,144.979,148.210,144.670,142.970,142.524,146.142,146.522,148.128,148.798,150.181,152.388,155.694,160.662,155.520,158.262,154.338,158.196,160.371,154.856,150.636,145.899,141.242,140.834,141.119,139.104,134.437,129.425,123.155,119.273,120.472,121.523,121.983,123.658,124.794,124.827,120.382,117.395,115.790,114.283,117.271,117.448,118.764,120.550,123.554,125.412,124.182,119.828,115.361,114.226,115.214,115.864,114.276,113.469,114.883,114.172,111.225,112.149,115.618,118.002,121.382,120.663,128.049) > 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 1.137138 69.202093 0.000000 0.000000 > m$fitted level slope sea Jan 1 183.04600 0.00000000 0.00000000 Feb 1 175.73240 -6.99649880 -0.01840302 Mar 1 175.80179 -0.05530959 -0.03379427 Apr 1 171.06263 -4.66446746 -0.03363156 May 1 170.49863 -0.62925201 -0.03363383 Jun 1 170.13563 -0.36723695 -0.03363383 Jul 1 156.42263 -13.50062473 -0.03363383 Aug 1 124.32463 -31.80205419 -0.03363383 Sep 1 99.39363 -25.04034122 -0.03363383 Oct 1 86.70863 -12.88161438 -0.03363383 Nov 1 85.08963 -1.79822548 -0.03363383 Dec 1 128.26963 42.46424752 -0.03363383 Jan 2 164.05002 35.90133511 0.20697543 Feb 2 162.56603 0.66617069 -0.16502635 Mar 2 152.92489 -9.46182068 -0.14588667 Apr 2 156.15126 3.02409587 -0.14626229 May 2 153.53326 -2.52821558 -0.14625963 Jun 2 153.33626 -0.23409736 -0.14625964 Jul 2 148.98626 -4.28450236 -0.14625964 Aug 2 144.35726 -4.62351790 -0.14625964 Sep 2 146.09926 1.64070354 -0.14625964 Oct 2 145.68826 -0.37835060 -0.14625964 Nov 2 150.41726 4.64772514 -0.14625964 Dec 2 147.63526 -2.66376848 -0.14625964 Jan 3 142.27208 -5.31513399 1.55192424 Feb 3 134.90214 -7.26218573 -0.14813974 Mar 3 131.89109 -3.08425308 -0.15509366 Apr 3 126.45903 -5.39463975 -0.15503245 May 3 125.66603 -0.86622732 -0.15503436 Jun 3 125.65003 -0.02952993 -0.15503437 Jul 3 130.28803 4.56372415 -0.15503437 Aug 3 126.41203 -3.74169605 -0.15503437 Sep 3 110.47803 -15.73998008 -0.15503437 Oct 3 98.57203 -11.96701132 -0.15503437 Nov 3 105.90403 7.02488884 -0.15503437 Dec 3 120.82003 14.79042630 -0.15503437 Jan 4 122.61114 2.01957373 1.46386087 Feb 4 127.36970 4.62580753 -0.12469596 Mar 4 146.87745 19.25354797 -0.14644931 Apr 4 145.12496 -1.41773363 -0.14595999 May 4 148.35596 3.15702326 -0.14596171 Jun 4 144.81596 -3.43342819 -0.14596167 Jul 4 143.11596 -1.72758458 -0.14596167 Aug 4 142.66996 -0.46639425 -0.14596167 Sep 4 146.28796 3.55300377 -0.14596167 Oct 4 146.66796 0.43049299 -0.14596167 Nov 4 148.27396 1.58729379 -0.14596167 Dec 4 148.94396 0.68459718 -0.14596167 Jan 5 148.59262 -0.33331721 1.58837941 Feb 5 152.52064 3.73508763 -0.13263588 Mar 5 155.82607 3.31272320 -0.13206821 Apr 5 160.79410 4.94165917 -0.13210305 May 5 155.65210 -4.98153560 -0.13209968 Jun 5 158.39410 2.61909298 -0.13209972 Jul 5 154.47010 -3.81987773 -0.13209972 Aug 5 158.32810 3.73581955 -0.13209972 Sep 5 160.50310 2.19983781 -0.13209972 Oct 5 154.98810 -5.39223139 -0.13209972 Nov 5 150.76810 -4.23865408 -0.13209972 Dec 5 146.03110 -4.72906967 -0.13209972 Jan 6 139.81127 -6.19410253 1.43072696 Feb 6 140.94564 0.82179960 -0.11164447 Mar 6 141.23000 0.29342244 -0.11099670 Apr 6 139.21495 -1.97826549 -0.11095237 May 6 134.54795 -4.62421333 -0.11095155 Jun 6 129.53595 -5.00582903 -0.11095154 Jul 6 123.26595 -6.24988285 -0.11095154 Aug 6 119.38395 -3.91968086 -0.11095154 Sep 6 120.58295 1.11754484 -0.11095154 Oct 6 121.63395 1.05205895 -0.11095154 Nov 6 122.09395 0.46942162 -0.11095154 Dec 6 123.76895 1.65581526 -0.11095154 Jan 7 123.59853 -0.13912075 1.19546718 Feb 7 124.93228 1.27415591 -0.10528006 Mar 7 120.48094 -4.35516252 -0.09893603 Apr 7 117.49396 -3.00877192 -0.09896018 May 7 115.88896 -1.62733866 -0.09896057 Jun 7 114.38196 -1.50891499 -0.09896057 Jul 7 117.36996 2.91643919 -0.09896057 Aug 7 117.54696 0.22059355 -0.09896057 Sep 7 118.86296 1.29856846 -0.09896057 Oct 7 120.64896 1.77824335 -0.09896057 Nov 7 123.65296 2.98449415 -0.09896057 Dec 7 125.51096 1.87592625 -0.09896057 Jan 8 123.14680 -2.29191681 1.03520404 Feb 8 119.92409 -3.18677370 -0.09609460 Mar 8 115.45578 -4.44687038 -0.09478058 Apr 8 114.32083 -1.18770266 -0.09483467 May 8 115.30884 0.95337737 -0.09483524 Jun 8 115.95884 0.65482774 -0.09483524 Jul 8 114.37084 -1.55230919 -0.09483524 Aug 8 113.56384 -0.81886034 -0.09483524 Sep 8 114.97784 1.37846780 -0.09483524 Oct 8 114.26684 -0.67774965 -0.09483524 Nov 8 111.31984 -2.91088871 -0.09483524 Dec 8 112.24384 0.86297422 -0.09483524 Jan 9 114.55791 2.28947404 1.06008701 Feb 9 118.09590 3.49178689 -0.09389638 Mar 9 121.47579 3.38175914 -0.09378960 Apr 9 120.75673 -0.65374353 -0.09372728 May 9 128.14273 7.25806106 -0.09372923 > m$resid Jan Feb Mar Apr May Jun 1 0.000000000 -0.859844815 0.823251470 -0.554064735 0.485073013 0.031496815 2 -0.799375597 -4.228955961 -1.203831829 1.500927056 -0.667443034 0.275775817 3 -0.322373128 -0.233764203 0.497271522 -0.277730755 0.544360193 0.100579344 4 -1.550689185 0.312977998 1.742899438 -2.484888025 0.549931257 -0.792237782 5 -0.123468221 0.488649179 -0.050368094 0.195813904 -1.192866661 0.913671116 6 -0.177549082 0.842773919 -0.063054969 -0.273078973 -0.318069235 -0.045874001 7 -0.217376298 0.169784061 -0.672184171 0.161849262 0.166062012 0.014235703 8 -0.504447654 -0.107511932 -0.150541024 0.391783750 0.257379105 -0.035888633 9 0.172566456 0.144459953 -0.013150469 -0.485106835 0.951077566 Jul Aug Sep Oct Nov Dec 1 -1.578763775 -2.200013764 0.812825120 1.461599841 1.332333531 5.320789287 2 -0.486898948 -0.040753038 0.753021692 -0.242710380 0.604184271 -0.878914220 3 0.552154804 -0.998394075 -1.442313021 0.453548354 2.283015205 0.933494803 4 0.205059360 0.151607615 0.483171600 -0.375356837 0.139058956 -0.108513107 5 -0.774028292 0.908269927 -0.184640275 -0.912642193 0.138671460 -0.058952829 6 -0.149547637 0.280113446 0.605524606 -0.007872055 -0.070038799 0.142616310 7 0.531971567 -0.324067447 0.129583301 0.057661691 0.145003338 -0.133260883 8 -0.265319800 0.088167843 0.264140687 -0.247177779 -0.268445516 0.453655843 9 > 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/wessaorg/rcomp/tmp/18k7y1324143155.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/wessaorg/rcomp/tmp/2bzys1324143155.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/wessaorg/rcomp/tmp/39cdo1324143155.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/wessaorg/rcomp/tmp/4l4461324143155.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/wessaorg/rcomp/tmp/5cz7d1324143155.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/wessaorg/rcomp/createtable file can be downloaded at http://www.wessa.net/cretab > load(file="/var/wessaorg/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/wessaorg/rcomp/tmp/69zbf1324143155.tab") > > try(system("convert tmp/18k7y1324143155.ps tmp/18k7y1324143155.png",intern=TRUE)) character(0) > try(system("convert tmp/2bzys1324143155.ps tmp/2bzys1324143155.png",intern=TRUE)) character(0) > try(system("convert tmp/39cdo1324143155.ps tmp/39cdo1324143155.png",intern=TRUE)) character(0) > try(system("convert tmp/4l4461324143155.ps tmp/4l4461324143155.png",intern=TRUE)) character(0) > try(system("convert tmp/5cz7d1324143155.ps tmp/5cz7d1324143155.png",intern=TRUE)) character(0) > > > proc.time() user system elapsed 2.053 0.383 2.449