R version 2.8.0 (2008-10-20) Copyright (C) 2008 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. > y <- c(116222,110924,103753,99983,93302,91496,119321,139261,133739,123913,113438,109416,109406,105645,101328,97686,93093,91382,122257,139183,139887,131822,116805,113706,113012,110452,107005,102841,98173,98181,137277,147579,146571,138920,130340,128140,127059,122860,117702,113537,108366,111078,150739,159129,157928,147768,137507,136919,136151,133001,125554,119647,114158,116193,152803,161761,160942,149470,139208,134588,130322,126611,122401,117352,112135,112879,148729,157230,157221,146681,136524,132111) > x <- c(492865,480961,461935,456608,441977,439148,488180,520564,501492,485025,464196,460170,467037,460070,447988,442867,436087,431328,484015,509673,512927,502831,470984,471067,476049,474605,470439,461251,454724,455626,516847,525192,522975,518585,509239,512238,519164,517009,509933,509127,500857,506971,569323,579714,577992,565464,547344,554788,562325,560854,555332,543599,536662,542722,593530,610763,612613,611324,594167,595454,590865,589379,584428,573100,567456,569028,620735,628884,628232,612117,595404,597141) > par2 = '36' > par1 = '0' > par1 <- as.numeric(par1) > par2 <- as.numeric(par2) > x <- as.ts(x) > y <- as.ts(y) > mylm <- lm(y~x) > cbind(mylm$resid) [,1] 1 -245.14170 2 -2277.98896 3 -4230.34011 4 -6539.19523 5 -9207.05261 6 -10237.08512 7 4138.90875 8 15196.28894 9 14905.55515 10 9596.29491 11 4834.48933 12 1916.78242 13 23.23035 14 -1826.78854 15 -2829.81211 16 -5067.17106 17 -7800.48225 18 -8206.13431 19 8217.32820 20 18105.58552 21 17917.04460 22 12621.28033 23 6339.60620 24 3217.84010 25 1157.32541 26 -1006.59928 27 -3310.90554 28 -4954.72548 29 -7832.43214 30 -8071.84207 31 14231.82630 32 22244.87304 33 21844.97482 34 15398.10961 35 9381.62754 36 6359.03068 37 3378.29548 38 -229.60874 39 -3446.72998 40 -7390.65193 41 -10293.27045 42 -9258.28186 43 13300.16408 44 18840.01269 45 18111.34074 46 11387.65067 47 6096.79244 48 3466.97482 49 631.64819 50 -2114.87065 51 -8047.23915 52 -10735.99007 53 -14322.23767 54 -13949.43740 55 8724.41674 56 12955.57054 57 11629.13333 58 510.69364 59 -5045.30623 60 -10018.31796 61 -13025.59938 62 -16329.00387 63 -19180.99218 64 -21122.83070 65 -24791.73577 66 -24478.92026 67 -2811.65317 68 3454.15449 69 3623.99182 70 -2495.81864 71 -8068.60344 72 -12958.04584 > library(lattice) > postscript(file="/var/www/html/freestat/rcomp/tmp/1py4w1257873570.ps",horizontal=F,pagecentre=F,paper="special",width=8.3333333333333,height=5.5555555555556) > plot(y,type='l',main='Run Sequence Plot of Y[t]',xlab='time or index',ylab='value') > grid() > dev.off() null device 1 > postscript(file="/var/www/html/freestat/rcomp/tmp/226gl1257873570.ps",horizontal=F,pagecentre=F,paper="special",width=8.3333333333333,height=5.5555555555556) > plot(x,type='l',main='Run Sequence Plot of X[t]',xlab='time or index',ylab='value') > grid() > dev.off() null device 1 > postscript(file="/var/www/html/freestat/rcomp/tmp/38pym1257873570.ps",horizontal=F,pagecentre=F,paper="special",width=8.3333333333333,height=5.5555555555556) > plot(x,y,main='Scatter Plot',xlab='X[t]',ylab='Y[t]') > grid() > dev.off() null device 1 > postscript(file="/var/www/html/freestat/rcomp/tmp/4ri5k1257873570.ps",horizontal=F,pagecentre=F,paper="special",width=8.3333333333333,height=5.5555555555556) > plot(mylm$resid,type='l',main='Run Sequence Plot of e[t]',xlab='time or index',ylab='value') > grid() > dev.off() null device 1 > postscript(file="/var/www/html/freestat/rcomp/tmp/5000m1257873570.ps",horizontal=F,pagecentre=F,paper="special",width=8.3333333333333,height=5.5555555555556) > hist(mylm$resid,main='Histogram of e[t]') > dev.off() null device 1 > postscript(file="/var/www/html/freestat/rcomp/tmp/6q44e1257873570.ps",horizontal=F,pagecentre=F,paper="special",width=8.3333333333333,height=5.5555555555556) > if (par1 > 0) + { + densityplot(~mylm$resid,col='black',main=paste('Density Plot of e[t] bw = ',par1),bw=par1) + } else { + densityplot(~mylm$resid,col='black',main='Density Plot of e[t]') + } > dev.off() null device 1 > postscript(file="/var/www/html/freestat/rcomp/tmp/7pbyi1257873570.ps",horizontal=F,pagecentre=F,paper="special",width=8.3333333333333,height=5.5555555555556) > qqnorm(mylm$resid,main='QQ plot of e[t]') > qqline(mylm$resid) > grid() > dev.off() null device 1 > if (par2 > 0) + { + postscript(file="/var/www/html/freestat/rcomp/tmp/8e3bu1257873570.ps",horizontal=F,pagecentre=F,paper="special",width=8.3333333333333,height=5.5555555555556) + acf(mylm$resid,lag.max=par2,main='Residual Autocorrelation Function') + grid() + dev.off() + } null device 1 > summary(x) Min. 1st Qu. Median Mean 3rd Qu. Max. 431300 473700 516900 524500 570300 628900 > > #Note: the /var/www/html/freestat/rcomp/createtable file can be downloaded at http://www.wessa.net/cretab > load(file="/var/www/html/freestat/rcomp/createtable") > > a<-table.start() > a<-table.row.start(a) > a<-table.element(a,'Model: Y[t] = c + b X[t] + e[t]',2,TRUE) > a<-table.row.end(a) > a<-table.row.start(a) > a<-table.element(a,'c',1,TRUE) > a<-table.element(a,mylm$coeff[[1]]) > a<-table.row.end(a) > a<-table.row.start(a) > a<-table.element(a,'b',1,TRUE) > a<-table.element(a,mylm$coeff[[2]]) > a<-table.row.end(a) > a<-table.end(a) > table.save(a,file="/var/www/html/freestat/rcomp/tmp/9o5t81257873570.tab") > a<-table.start() > a<-table.row.start(a) > a<-table.element(a,'Descriptive Statistics about e[t]',2,TRUE) > a<-table.row.end(a) > a<-table.row.start(a) > a<-table.element(a,'# observations',header=TRUE) > a<-table.element(a,length(mylm$resid)) > a<-table.row.end(a) > a<-table.row.start(a) > a<-table.element(a,'minimum',header=TRUE) > a<-table.element(a,min(mylm$resid)) > a<-table.row.end(a) > a<-table.row.start(a) > a<-table.element(a,'Q1',header=TRUE) > a<-table.element(a,quantile(mylm$resid,0.25)) > a<-table.row.end(a) > a<-table.row.start(a) > a<-table.element(a,'median',header=TRUE) > a<-table.element(a,median(mylm$resid)) > a<-table.row.end(a) > a<-table.row.start(a) > a<-table.element(a,'mean',header=TRUE) > a<-table.element(a,mean(mylm$resid)) > a<-table.row.end(a) > a<-table.row.start(a) > a<-table.element(a,'Q3',header=TRUE) > a<-table.element(a,quantile(mylm$resid,0.75)) > a<-table.row.end(a) > a<-table.row.start(a) > a<-table.element(a,'maximum',header=TRUE) > a<-table.element(a,max(mylm$resid)) > a<-table.row.end(a) > a<-table.end(a) > table.save(a,file="/var/www/html/freestat/rcomp/tmp/10pnw51257873570.tab") > > system("convert tmp/1py4w1257873570.ps tmp/1py4w1257873570.png") > system("convert tmp/226gl1257873570.ps tmp/226gl1257873570.png") > system("convert tmp/38pym1257873570.ps tmp/38pym1257873570.png") > system("convert tmp/4ri5k1257873570.ps tmp/4ri5k1257873570.png") > system("convert tmp/5000m1257873570.ps tmp/5000m1257873570.png") > system("convert tmp/6q44e1257873570.ps tmp/6q44e1257873570.png") > system("convert tmp/7pbyi1257873570.ps tmp/7pbyi1257873570.png") > system("convert tmp/8e3bu1257873570.ps tmp/8e3bu1257873570.png") > > > proc.time() user system elapsed 2.651 2.030 3.082