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(8.4,8.4,8.4,8.6,8.9,8.8,8.3,7.5,7.2,7.4,8.8,9.3,9.3,8.7,8.2,8.3,8.5,8.6,8.5,8.2,8.1,7.9,8.6,8.7,8.7,8.5,8.4,8.5,8.7,8.7,8.6,8.5,8.3,8,8.2,8.1,8.1,8,7.9,7.9,8,8,7.9,8,7.7,7.2,7.5,7.3,7,7,7,7.2,7.3,7.1,6.8,6.4,6.1,6.5,7.7,7.9,7.5,6.9,6.6,6.9,7.7,8,8,7.7,7.3,7.4,8.1,8.3,8.2) > 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.0000000000 0.1179061165 0.0009824949 0.0000000000 > m$fitted level slope sea Jan 1 8.400000 0.000000000 0.0000000000 Feb 1 8.400000 0.000000000 0.0000000000 Mar 1 8.400000 0.000000000 0.0000000000 Apr 1 8.594503 0.186336498 0.0054971182 May 1 8.898188 0.298820879 0.0018124251 Jun 1 8.812189 -0.070103690 -0.0121887942 Jul 1 8.309828 -0.484497369 -0.0098282071 Aug 1 7.506194 -0.790443779 -0.0061943703 Sep 1 7.184198 -0.341356654 0.0158018200 Oct 1 7.387840 0.181117624 0.0121600815 Nov 1 8.769403 1.331950182 0.0305968667 Dec 1 9.331418 0.593835381 -0.0314179686 Jan 2 9.312060 0.006142299 -0.0120602462 Feb 2 8.714171 -0.573184008 -0.0141712139 Mar 2 8.192914 -0.524590012 0.0070862224 Apr 2 8.278707 0.046084776 0.0212934814 May 2 8.490827 0.201233414 0.0091734296 Jun 2 8.610260 0.124775166 -0.0102601131 Jul 2 8.499062 -0.095789375 0.0009383713 Aug 2 8.220600 -0.266530263 -0.0205995037 Sep 2 8.067124 -0.160859902 0.0328757339 Oct 2 7.951918 -0.118188339 -0.0519178929 Nov 2 8.519741 0.523020367 0.0802590403 Dec 2 8.742425 0.242306274 -0.0424253381 Jan 3 8.697978 -0.025641721 0.0020215728 Feb 3 8.501323 -0.185569563 -0.0013225125 Mar 3 8.418115 -0.091154706 -0.0181146935 Apr 3 8.483651 0.053599204 0.0163488407 May 3 8.687386 0.192128801 0.0126138218 Jun 3 8.711844 0.037367502 -0.0118439488 Jul 3 8.593042 -0.106783564 0.0069584733 Aug 3 8.521458 -0.074293406 -0.0214578087 Sep 3 8.255308 -0.251378878 0.0446918241 Oct 3 8.099858 -0.162835389 -0.0998575888 Nov 3 8.101773 -0.010765202 0.0982265874 Dec 3 8.129558 0.024815043 -0.0295580938 Jan 4 8.088101 -0.036339453 0.0118985801 Feb 4 8.004097 -0.080343798 -0.0040971136 Mar 4 7.924103 -0.080023604 -0.0241028133 Apr 4 7.898612 -0.029968925 0.0013879897 May 4 7.973795 0.066426835 0.0262052415 Jun 4 8.002771 0.032085476 -0.0027708730 Jul 4 7.910181 -0.082247699 -0.0101805591 Aug 4 7.993013 0.069135880 0.0069869848 Sep 4 7.677261 -0.283812957 0.0227386764 Oct 4 7.314196 -0.356489541 -0.1141959239 Nov 4 7.382784 0.033320747 0.1172158796 Dec 4 7.332400 -0.043431179 -0.0323998107 Jan 5 7.001452 -0.307046240 -0.0014515762 Feb 5 6.984670 -0.040900703 0.0153302726 Mar 5 7.020540 0.029174394 -0.0205397243 Apr 5 7.201629 0.168063169 -0.0016294266 May 5 7.278294 0.084587780 0.0217058220 Jun 5 7.104478 -0.151449040 -0.0044776687 Jul 5 6.844911 -0.250218854 -0.0449109586 Aug 5 6.370486 -0.455039515 0.0295141464 Sep 5 6.046628 -0.335201367 0.0533715829 Oct 5 6.606302 0.482295506 -0.1063015882 Nov 5 7.561006 0.913854115 0.1389940990 Dec 5 7.931980 0.417984604 -0.0319803309 Jan 6 7.566743 -0.297445557 -0.0667427444 Feb 6 6.903715 -0.631307124 -0.0037148709 Mar 6 6.614648 -0.319647775 -0.0146481636 Apr 6 6.868996 0.203719305 0.0310039867 May 6 7.635114 0.716126230 0.0648855123 Jun 6 8.026661 0.420394583 -0.0266609509 Jul 6 8.037309 0.047016797 -0.0373088679 Aug 6 7.663819 -0.336168985 0.0361806245 Sep 6 7.307182 -0.354820773 -0.0071818321 Oct 6 7.539107 0.179848045 -0.1391067898 Nov 6 7.935712 0.377362761 0.1642879470 Dec 6 8.280997 0.348136769 0.0190028805 Jan 7 8.241435 -0.005136546 -0.0414345341 > m$resid Jan Feb Mar Apr May Jun 1 0.000000000 0.000000000 0.000000000 0.542139809 0.327604743 -1.074407860 2 -1.712190020 -1.694238640 0.142937770 1.660002719 0.451892591 -0.222666517 3 -0.780709101 -0.466614957 0.276352870 0.421298119 0.403448764 -0.450712780 4 -0.178202264 -0.128225093 0.000934917 0.145753331 0.280701723 -0.100014939 5 -0.768172856 0.775123402 0.204371834 0.404538840 -0.243056861 -0.687434973 6 -2.084625947 -0.972161148 0.908425294 1.524585347 1.491935988 -0.861286444 7 -1.029278638 Jul Aug Sep Oct Nov Dec 1 -1.206826938 -0.890998953 1.307863551 1.521586854 3.351536643 -2.149590560 2 -0.642344402 -0.497243774 0.307740764 0.124271170 1.867374453 -0.817517190 3 -0.419806818 0.094620183 -0.515720954 0.257862671 0.442870738 0.103619749 4 -0.332968049 0.440870368 -1.027882786 -0.211653986 1.135237907 -0.223525682 5 -0.287643357 -0.596493447 0.349001216 2.380775477 1.256818997 -1.444159268 6 -1.087376664 -1.115939674 -0.054319087 1.557104625 0.575216160 -0.085119325 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/www/html/rcomp/tmp/146uw1259697275.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/24tpx1259697275.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/344i21259697275.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/4erhh1259697275.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/5k4v41259697275.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/6u45l1259697275.tab") > system("convert tmp/146uw1259697275.ps tmp/146uw1259697275.png") > system("convert tmp/24tpx1259697275.ps tmp/24tpx1259697275.png") > system("convert tmp/344i21259697275.ps tmp/344i21259697275.png") > system("convert tmp/4erhh1259697275.ps tmp/4erhh1259697275.png") > system("convert tmp/5k4v41259697275.ps tmp/5k4v41259697275.png") > > > proc.time() user system elapsed 1.450 0.821 1.758