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(103.63,103.64,103.66,103.77,103.88,103.91,103.91,103.92,104.05,104.23,104.30,104.31,104.31,104.34,104.55,104.65,104.73,104.75,104.75,104.76,104.94,105.29,105.38,105.43,105.43,105.42,105.52,105.69,105.72,105.74,105.74,105.74,105.95,106.17,106.34,106.37,106.37,106.36,106.44,106.29,106.23,106.23,106.23,106.23,106.34,106.44,106.44,106.48,106.50,106.57,106.40,106.37,106.25,106.21,106.21,106.24,106.19,106.08,106.13,106.09) > 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 0.0065962768 0.0002806596 0.0000000000 0.0000000000 > m$fitted level slope sea Jan 1 103.6300 0.0000000000 0.0000000000 Feb 1 103.6395 0.0006489105 0.0005129289 Mar 1 103.6595 0.0021528282 0.0005473931 Apr 1 103.7693 0.0137307163 0.0007188477 May 1 103.8791 0.0262805009 0.0008519442 Jun 1 103.9091 0.0268288267 0.0008563285 Jul 1 103.9092 0.0225458493 0.0008297530 Aug 1 103.9192 0.0204358372 0.0008194157 Sep 1 104.0491 0.0395058701 0.0008939794 Oct 1 104.2290 0.0645192108 0.0009725696 Nov 1 104.2990 0.0655097086 0.0009750814 Dec 1 104.3090 0.0553781357 0.0009542850 Jan 2 104.3336 0.0500146811 -0.0235569550 Feb 2 104.3389 0.0419403927 0.0011294292 Mar 2 104.5490 0.0730644552 0.0010368585 Apr 2 104.6490 0.0780609571 0.0010247738 May 2 104.7290 0.0784210348 0.0010240654 Jun 2 104.7490 0.0675645707 0.0010414425 Jul 2 104.7489 0.0550030057 0.0010578030 Aug 2 104.7589 0.0466334496 0.0010666736 Sep 2 104.9390 0.0714418079 0.0010452755 Oct 2 105.2890 0.1232653148 0.0010088968 Nov 2 105.3790 0.1170760069 0.0010124329 Dec 2 105.4290 0.1045951654 0.0010182362 Jan 3 105.4567 0.0904927124 -0.0266621681 Feb 3 105.4194 0.0672061618 0.0005716906 Mar 3 105.5194 0.0733139288 0.0005569784 Apr 3 105.6895 0.0913163920 0.0005216725 May 3 105.7195 0.0799017147 0.0005398946 Jun 3 105.7394 0.0687517715 0.0005543828 Jul 3 105.7394 0.0559555497 0.0005679165 Aug 3 105.7394 0.0455415516 0.0005768813 Sep 3 105.9494 0.0761481399 0.0005554365 Oct 3 106.1695 0.1029190961 0.0005401696 Nov 3 106.3395 0.1154027140 0.0005343752 Dec 3 106.3695 0.0995096110 0.0005403793 Jan 4 106.3868 0.0843515477 -0.0167748524 Feb 4 106.3597 0.0639284105 0.0003326706 Mar 4 106.4397 0.0669212106 0.0003272403 Apr 4 106.2896 0.0265359572 0.0003869241 May 4 106.2296 0.0104275588 0.0004063016 Jun 4 106.2296 0.0084866938 0.0004082020 Jul 4 106.2296 0.0069071828 0.0004094608 Aug 4 106.2296 0.0056217002 0.0004102947 Sep 4 106.3396 0.0250467852 0.0004000389 Oct 4 106.4396 0.0389955155 0.0003940449 Nov 4 106.4396 0.0317385746 0.0003965831 Dec 4 106.4796 0.0332759866 0.0003961454 Jan 5 106.5051 0.0318411940 -0.0051180723 Feb 5 106.5693 0.0377789513 0.0007391029 Mar 5 106.3992 -0.0009077274 0.0007949790 Apr 5 106.3692 -0.0063234861 0.0008013458 May 5 106.2492 -0.0274827513 0.0008215928 Jun 5 106.2092 -0.0298124749 0.0008234073 Jul 5 106.2092 -0.0242640362 0.0008198900 Aug 5 106.2392 -0.0141652357 0.0008146793 Sep 5 106.1892 -0.0208341060 0.0008174799 Oct 5 106.0792 -0.0374276776 0.0008231518 Nov 5 106.1292 -0.0211577482 0.0008186255 Dec 5 106.0892 -0.0246641900 0.0008194194 > m$resid Jan Feb Mar Apr May Jun 1 0.00000000 0.06960653 0.22881683 1.25458131 1.10537018 0.04228660 2 -0.40240433 -0.42032620 1.86784282 0.29931073 0.02154421 -0.64904824 3 -0.92733967 -1.28507599 0.36424112 1.07392673 -0.68107760 -0.66537333 4 -0.96514502 -1.15357189 0.17851045 -2.40945135 -0.96121421 -0.11582736 5 -0.08991660 0.33959071 -2.30788001 -0.32314561 -1.26269160 -0.13903987 Jul Aug Sep Oct Nov Dec 1 -0.30282010 -0.14088483 1.22601287 1.56831943 0.06108146 -0.61795715 2 -0.75059088 -0.49993119 1.48151251 3.09433630 -0.36952063 -0.74509433 3 -0.76368471 -0.62154998 1.82679941 1.59790640 0.74513600 -0.94865615 4 -0.09426930 -0.07672464 1.15943280 0.83258115 -0.43316325 0.09176823 5 0.33115484 0.60276342 -0.39805262 -0.99045662 0.97114959 -0.20930053 > 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/html/rcomp/tmp/11au41259502258.ps",horizontal=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/290bd1259502258.ps",horizontal=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/39n3o1259502258.ps",horizontal=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/4v7so1259502258.ps",horizontal=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/html/rcomp/tmp/5vrf61259502258.ps",horizontal=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/692x51259502258.tab") > system("convert tmp/11au41259502258.ps tmp/11au41259502258.png") > system("convert tmp/290bd1259502258.ps tmp/290bd1259502258.png") > system("convert tmp/39n3o1259502258.ps tmp/39n3o1259502258.png") > system("convert tmp/4v7so1259502258.ps tmp/4v7so1259502258.png") > system("convert tmp/5vrf61259502258.ps tmp/5vrf61259502258.png") > > > proc.time() user system elapsed 1.329 0.801 1.521