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 <- array(list(105.7,0,105.7,0,111.1,0,82.4,0,60,0,107.3,0,99.3,0,113.5,0,108.9,0,100.2,0,103.9,0,138.7,0,120.2,0,100.2,0,143.2,0,70.9,0,85.2,0,133,0,136.6,0,117.9,0,106.3,0,122.3,0,125.5,0,148.4,0,126.3,0,99.6,0,140.4,0,80.3,0,92.6,0,138.5,0,110.9,0,119.6,0,105,0,109,0,129.4,0,148.6,0,101.4,0,134.8,0,143.7,0,81.6,0,90.3,0,141.5,0,140.7,0,140.2,0,100.2,0,125.7,0,119.6,0,134.7,0,109,0,116.3,0,146.9,0,97.4,0,89.4,0,132.1,0,139.8,0,129,0,112.5,0,121.9,0,121.7,0,123.1,0,131.6,0,119.3,0,132.5,0,98.3,0,85.1,0,131.7,0,129.3,0,90.7,1,78.6,1,68.9,1,79.1,1,83.5,1,74.1,1,59.7,1,93.3,1,61.3,1,56.6,1),dim=c(2,77),dimnames=list(c('Y','X'),1:77))
> y <- array(NA,dim=c(2,77),dimnames=list(c('Y','X'),1:77))
> for (i in 1:dim(x)[1])
+ {
+ for (j in 1:dim(x)[2])
+ {
+ y[i,j] <- as.numeric(x[i,j])
+ }
+ }
> par3 = 'Linear Trend'
> par2 = 'Include Monthly Dummies'
> par1 = '1'
> #'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!)
> library(lattice)
> library(lmtest)
Loading required package: zoo
Attaching package: 'zoo'
The following object(s) are masked from package:base :
as.Date.numeric
> n25 <- 25 #minimum number of obs. for Goldfeld-Quandt test
> par1 <- as.numeric(par1)
> x <- t(y)
> k <- length(x[1,])
> n <- length(x[,1])
> x1 <- cbind(x[,par1], x[,1:k!=par1])
> mycolnames <- c(colnames(x)[par1], colnames(x)[1:k!=par1])
> colnames(x1) <- mycolnames #colnames(x)[par1]
> x <- x1
> if (par3 == 'First Differences'){
+ x2 <- array(0, dim=c(n-1,k), dimnames=list(1:(n-1), paste('(1-B)',colnames(x),sep='')))
+ for (i in 1:n-1) {
+ for (j in 1:k) {
+ x2[i,j] <- x[i+1,j] - x[i,j]
+ }
+ }
+ x <- x2
+ }
> if (par2 == 'Include Monthly Dummies'){
+ x2 <- array(0, dim=c(n,11), dimnames=list(1:n, paste('M', seq(1:11), sep ='')))
+ for (i in 1:11){
+ x2[seq(i,n,12),i] <- 1
+ }
+ x <- cbind(x, x2)
+ }
> if (par2 == 'Include Quarterly Dummies'){
+ x2 <- array(0, dim=c(n,3), dimnames=list(1:n, paste('Q', seq(1:3), sep ='')))
+ for (i in 1:3){
+ x2[seq(i,n,4),i] <- 1
+ }
+ x <- cbind(x, x2)
+ }
> k <- length(x[1,])
> if (par3 == 'Linear Trend'){
+ x <- cbind(x, c(1:n))
+ colnames(x)[k+1] <- 't'
+ }
> x
Y X M1 M2 M3 M4 M5 M6 M7 M8 M9 M10 M11 t
1 105.7 0 1 0 0 0 0 0 0 0 0 0 0 1
2 105.7 0 0 1 0 0 0 0 0 0 0 0 0 2
3 111.1 0 0 0 1 0 0 0 0 0 0 0 0 3
4 82.4 0 0 0 0 1 0 0 0 0 0 0 0 4
5 60.0 0 0 0 0 0 1 0 0 0 0 0 0 5
6 107.3 0 0 0 0 0 0 1 0 0 0 0 0 6
7 99.3 0 0 0 0 0 0 0 1 0 0 0 0 7
8 113.5 0 0 0 0 0 0 0 0 1 0 0 0 8
9 108.9 0 0 0 0 0 0 0 0 0 1 0 0 9
10 100.2 0 0 0 0 0 0 0 0 0 0 1 0 10
11 103.9 0 0 0 0 0 0 0 0 0 0 0 1 11
12 138.7 0 0 0 0 0 0 0 0 0 0 0 0 12
13 120.2 0 1 0 0 0 0 0 0 0 0 0 0 13
14 100.2 0 0 1 0 0 0 0 0 0 0 0 0 14
15 143.2 0 0 0 1 0 0 0 0 0 0 0 0 15
16 70.9 0 0 0 0 1 0 0 0 0 0 0 0 16
17 85.2 0 0 0 0 0 1 0 0 0 0 0 0 17
18 133.0 0 0 0 0 0 0 1 0 0 0 0 0 18
19 136.6 0 0 0 0 0 0 0 1 0 0 0 0 19
20 117.9 0 0 0 0 0 0 0 0 1 0 0 0 20
21 106.3 0 0 0 0 0 0 0 0 0 1 0 0 21
22 122.3 0 0 0 0 0 0 0 0 0 0 1 0 22
23 125.5 0 0 0 0 0 0 0 0 0 0 0 1 23
24 148.4 0 0 0 0 0 0 0 0 0 0 0 0 24
25 126.3 0 1 0 0 0 0 0 0 0 0 0 0 25
26 99.6 0 0 1 0 0 0 0 0 0 0 0 0 26
27 140.4 0 0 0 1 0 0 0 0 0 0 0 0 27
28 80.3 0 0 0 0 1 0 0 0 0 0 0 0 28
29 92.6 0 0 0 0 0 1 0 0 0 0 0 0 29
30 138.5 0 0 0 0 0 0 1 0 0 0 0 0 30
31 110.9 0 0 0 0 0 0 0 1 0 0 0 0 31
32 119.6 0 0 0 0 0 0 0 0 1 0 0 0 32
33 105.0 0 0 0 0 0 0 0 0 0 1 0 0 33
34 109.0 0 0 0 0 0 0 0 0 0 0 1 0 34
35 129.4 0 0 0 0 0 0 0 0 0 0 0 1 35
36 148.6 0 0 0 0 0 0 0 0 0 0 0 0 36
37 101.4 0 1 0 0 0 0 0 0 0 0 0 0 37
38 134.8 0 0 1 0 0 0 0 0 0 0 0 0 38
39 143.7 0 0 0 1 0 0 0 0 0 0 0 0 39
40 81.6 0 0 0 0 1 0 0 0 0 0 0 0 40
41 90.3 0 0 0 0 0 1 0 0 0 0 0 0 41
42 141.5 0 0 0 0 0 0 1 0 0 0 0 0 42
43 140.7 0 0 0 0 0 0 0 1 0 0 0 0 43
44 140.2 0 0 0 0 0 0 0 0 1 0 0 0 44
45 100.2 0 0 0 0 0 0 0 0 0 1 0 0 45
46 125.7 0 0 0 0 0 0 0 0 0 0 1 0 46
47 119.6 0 0 0 0 0 0 0 0 0 0 0 1 47
48 134.7 0 0 0 0 0 0 0 0 0 0 0 0 48
49 109.0 0 1 0 0 0 0 0 0 0 0 0 0 49
50 116.3 0 0 1 0 0 0 0 0 0 0 0 0 50
51 146.9 0 0 0 1 0 0 0 0 0 0 0 0 51
52 97.4 0 0 0 0 1 0 0 0 0 0 0 0 52
53 89.4 0 0 0 0 0 1 0 0 0 0 0 0 53
54 132.1 0 0 0 0 0 0 1 0 0 0 0 0 54
55 139.8 0 0 0 0 0 0 0 1 0 0 0 0 55
56 129.0 0 0 0 0 0 0 0 0 1 0 0 0 56
57 112.5 0 0 0 0 0 0 0 0 0 1 0 0 57
58 121.9 0 0 0 0 0 0 0 0 0 0 1 0 58
59 121.7 0 0 0 0 0 0 0 0 0 0 0 1 59
60 123.1 0 0 0 0 0 0 0 0 0 0 0 0 60
61 131.6 0 1 0 0 0 0 0 0 0 0 0 0 61
62 119.3 0 0 1 0 0 0 0 0 0 0 0 0 62
63 132.5 0 0 0 1 0 0 0 0 0 0 0 0 63
64 98.3 0 0 0 0 1 0 0 0 0 0 0 0 64
65 85.1 0 0 0 0 0 1 0 0 0 0 0 0 65
66 131.7 0 0 0 0 0 0 1 0 0 0 0 0 66
67 129.3 0 0 0 0 0 0 0 1 0 0 0 0 67
68 90.7 1 0 0 0 0 0 0 0 1 0 0 0 68
69 78.6 1 0 0 0 0 0 0 0 0 1 0 0 69
70 68.9 1 0 0 0 0 0 0 0 0 0 1 0 70
71 79.1 1 0 0 0 0 0 0 0 0 0 0 1 71
72 83.5 1 0 0 0 0 0 0 0 0 0 0 0 72
73 74.1 1 1 0 0 0 0 0 0 0 0 0 0 73
74 59.7 1 0 1 0 0 0 0 0 0 0 0 0 74
75 93.3 1 0 0 1 0 0 0 0 0 0 0 0 75
76 61.3 1 0 0 0 1 0 0 0 0 0 0 0 76
77 56.6 1 0 0 0 0 1 0 0 0 0 0 0 77
> k <- length(x[1,])
> df <- as.data.frame(x)
> (mylm <- lm(df))
Call:
lm(formula = df)
Coefficients:
(Intercept) X M1 M2 M3 M4
126.3838 -49.8979 -19.5699 -24.5135 0.2857 -48.4008
M5 M6 M7 M8 M9 M10
-50.5301 -5.4998 -10.3553 -9.9279 -26.7667 -20.9556
M11 t
-16.0278 0.2722
> (mysum <- summary(mylm))
Call:
lm(formula = df)
Residuals:
Min 1Q Median 3Q Max
-19.6160 -5.6846 -0.6846 8.8003 22.5860
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 126.38376 4.92436 25.665 < 2e-16 ***
X -49.89786 4.44272 -11.231 < 2e-16 ***
M1 -19.56988 5.84107 -3.350 0.001367 **
M2 -24.51352 5.83837 -4.199 8.59e-05 ***
M3 0.28571 5.83643 0.049 0.961112
M4 -48.40078 5.83524 -8.295 1.09e-11 ***
M5 -50.53013 5.83481 -8.660 2.52e-12 ***
M6 -5.49976 6.08428 -0.904 0.369477
M7 -10.35529 6.08501 -1.702 0.093729 .
M8 -9.92785 6.06024 -1.638 0.106367
M9 -26.76672 6.05769 -4.419 3.99e-05 ***
M10 -20.95559 6.05587 -3.460 0.000974 ***
M11 -16.02780 6.05478 -2.647 0.010242 *
t 0.27220 0.06643 4.097 0.000122 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 10.49 on 63 degrees of freedom
Multiple R-squared: 0.8427, Adjusted R-squared: 0.8103
F-statistic: 25.97 on 13 and 63 DF, p-value: < 2.2e-16
> if (n > n25) {
+ kp3 <- k + 3
+ nmkm3 <- n - k - 3
+ gqarr <- array(NA, dim=c(nmkm3-kp3+1,3))
+ numgqtests <- 0
+ numsignificant1 <- 0
+ numsignificant5 <- 0
+ numsignificant10 <- 0
+ for (mypoint in kp3:nmkm3) {
+ j <- 0
+ numgqtests <- numgqtests + 1
+ for (myalt in c('greater', 'two.sided', 'less')) {
+ j <- j + 1
+ gqarr[mypoint-kp3+1,j] <- gqtest(mylm, point=mypoint, alternative=myalt)$p.value
+ }
+ if (gqarr[mypoint-kp3+1,2] < 0.01) numsignificant1 <- numsignificant1 + 1
+ if (gqarr[mypoint-kp3+1,2] < 0.05) numsignificant5 <- numsignificant5 + 1
+ if (gqarr[mypoint-kp3+1,2] < 0.10) numsignificant10 <- numsignificant10 + 1
+ }
+ gqarr
+ }
[,1] [,2] [,3]
[1,] 0.9362056 0.12758887 0.06379443
[2,] 0.9092769 0.18144613 0.09072306
[3,] 0.9298979 0.14020423 0.07010212
[4,] 0.9081727 0.18365452 0.09182726
[5,] 0.8984991 0.20300171 0.10150085
[6,] 0.8571725 0.28565505 0.14282752
[7,] 0.8016936 0.39661286 0.19830643
[8,] 0.7687374 0.46252519 0.23126259
[9,] 0.7394637 0.52107262 0.26053631
[10,] 0.8316569 0.33668621 0.16834310
[11,] 0.7705449 0.45891013 0.22945506
[12,] 0.7637951 0.47240982 0.23620491
[13,] 0.7009595 0.59808103 0.29904051
[14,] 0.6300149 0.73997025 0.36998512
[15,] 0.8015869 0.39682617 0.19841308
[16,] 0.8046716 0.39065675 0.19532838
[17,] 0.8017146 0.39657072 0.19828536
[18,] 0.8132392 0.37352169 0.18676085
[19,] 0.7631690 0.47366207 0.23683103
[20,] 0.7908435 0.41831302 0.20915651
[21,] 0.9318287 0.13634266 0.06817133
[22,] 0.9757168 0.04856645 0.02428322
[23,] 0.9633190 0.07336192 0.03668096
[24,] 0.9829873 0.03402534 0.01701267
[25,] 0.9736708 0.05265843 0.02632921
[26,] 0.9623980 0.07520402 0.03760201
[27,] 0.9515659 0.09686815 0.04843408
[28,] 0.9424033 0.11519341 0.05759670
[29,] 0.9756129 0.04877427 0.02438714
[30,] 0.9659403 0.06811947 0.03405974
[31,] 0.9529764 0.09404711 0.04702356
[32,] 0.9581012 0.08379769 0.04189884
[33,] 0.9823255 0.03534897 0.01767449
[34,] 0.9692447 0.06151065 0.03075533
[35,] 0.9672778 0.06544431 0.03272216
[36,] 0.9493107 0.10137865 0.05068933
[37,] 0.9270952 0.14580964 0.07290482
[38,] 0.8987485 0.20250299 0.10125149
[39,] 0.8378819 0.32423617 0.16211808
[40,] 0.7636100 0.47277996 0.23638998
[41,] 0.7116459 0.57670821 0.28835411
[42,] 0.6342500 0.73149996 0.36574998
[43,] 0.4898338 0.97966764 0.51016618
[44,] 0.3963427 0.79268536 0.60365732
> postscript(file="/var/www/html/rcomp/tmp/1b32o1258731436.ps",horizontal=F,pagecentre=F,paper="special",width=8.3333333333333,height=5.5555555555556)
> plot(x[,1], type='l', main='Actuals and Interpolation', ylab='value of Actuals and Interpolation (dots)', xlab='time or index')
> points(x[,1]-mysum$resid)
> grid()
> dev.off()
null device
1
> postscript(file="/var/www/html/rcomp/tmp/2fjcb1258731436.ps",horizontal=F,pagecentre=F,paper="special",width=8.3333333333333,height=5.5555555555556)
> plot(mysum$resid, type='b', pch=19, main='Residuals', ylab='value of Residuals', xlab='time or index')
> grid()
> dev.off()
null device
1
> postscript(file="/var/www/html/rcomp/tmp/3hjpd1258731436.ps",horizontal=F,pagecentre=F,paper="special",width=8.3333333333333,height=5.5555555555556)
> hist(mysum$resid, main='Residual Histogram', xlab='values of Residuals')
> grid()
> dev.off()
null device
1
> postscript(file="/var/www/html/rcomp/tmp/4006s1258731436.ps",horizontal=F,pagecentre=F,paper="special",width=8.3333333333333,height=5.5555555555556)
> densityplot(~mysum$resid,col='black',main='Residual Density Plot', xlab='values of Residuals')
> dev.off()
null device
1
> postscript(file="/var/www/html/rcomp/tmp/5d0zx1258731436.ps",horizontal=F,pagecentre=F,paper="special",width=8.3333333333333,height=5.5555555555556)
> qqnorm(mysum$resid, main='Residual Normal Q-Q Plot')
> qqline(mysum$resid)
> grid()
> dev.off()
null device
1
> (myerror <- as.ts(mysum$resid))
Time Series:
Start = 1
End = 77
Frequency = 1
1 2 3 4 5 6
-1.3860806 3.2853480 -16.3860806 3.3282051 -17.2146520 -15.2172263
7 8 9 10 11 12
-18.6338930 -5.1335368 6.8331298 -7.9502035 -9.4502035 9.0497965
13 14 15 16 17 18
9.8474766 -5.4810948 12.4474766 -11.4382377 4.7189052 7.2163309
19 20 21 22 23 24
15.3996642 -3.9999796 0.9666870 10.8833537 8.8833537 15.4833537
25 26 27 28 29 30
12.6810338 -9.3475376 6.3810338 -5.3046805 8.8524624 9.4498881
31 32 33 34 35 36
-13.5667786 -5.5664225 -3.5997558 -5.6830891 9.5169109 12.4169109
37 38 39 40 41 42
-15.4854090 22.5860195 6.4145910 -7.2711233 3.2860195 9.1834453
43 44 45 46 47 48
12.9667786 11.7671347 -11.6661986 7.7504681 -3.5495319 -4.7495319
49 50 51 52 53 54
-11.1518519 0.8195767 6.3481481 5.2624339 -0.8804233 -3.4829976
55 56 57 58 59 60
8.8003358 -2.6993081 -2.6326414 0.6840252 -4.7159748 -19.6159748
61 62 63 64 65 66
8.1817053 0.5531339 -11.3182947 2.8959910 -8.4468661 -7.1494404
67 68 69 70 71 72
-4.9661070 5.6321123 10.0987790 -5.6845543 -0.6845543 -12.5845543
73 74 75 76 77
-2.6868742 -12.4154457 -3.8868742 12.5274115 9.6845543
> postscript(file="/var/www/html/rcomp/tmp/6qvdx1258731436.ps",horizontal=F,pagecentre=F,paper="special",width=8.3333333333333,height=5.5555555555556)
> dum <- cbind(lag(myerror,k=1),myerror)
> dum
Time Series:
Start = 0
End = 77
Frequency = 1
lag(myerror, k = 1) myerror
0 -1.3860806 NA
1 3.2853480 -1.3860806
2 -16.3860806 3.2853480
3 3.3282051 -16.3860806
4 -17.2146520 3.3282051
5 -15.2172263 -17.2146520
6 -18.6338930 -15.2172263
7 -5.1335368 -18.6338930
8 6.8331298 -5.1335368
9 -7.9502035 6.8331298
10 -9.4502035 -7.9502035
11 9.0497965 -9.4502035
12 9.8474766 9.0497965
13 -5.4810948 9.8474766
14 12.4474766 -5.4810948
15 -11.4382377 12.4474766
16 4.7189052 -11.4382377
17 7.2163309 4.7189052
18 15.3996642 7.2163309
19 -3.9999796 15.3996642
20 0.9666870 -3.9999796
21 10.8833537 0.9666870
22 8.8833537 10.8833537
23 15.4833537 8.8833537
24 12.6810338 15.4833537
25 -9.3475376 12.6810338
26 6.3810338 -9.3475376
27 -5.3046805 6.3810338
28 8.8524624 -5.3046805
29 9.4498881 8.8524624
30 -13.5667786 9.4498881
31 -5.5664225 -13.5667786
32 -3.5997558 -5.5664225
33 -5.6830891 -3.5997558
34 9.5169109 -5.6830891
35 12.4169109 9.5169109
36 -15.4854090 12.4169109
37 22.5860195 -15.4854090
38 6.4145910 22.5860195
39 -7.2711233 6.4145910
40 3.2860195 -7.2711233
41 9.1834453 3.2860195
42 12.9667786 9.1834453
43 11.7671347 12.9667786
44 -11.6661986 11.7671347
45 7.7504681 -11.6661986
46 -3.5495319 7.7504681
47 -4.7495319 -3.5495319
48 -11.1518519 -4.7495319
49 0.8195767 -11.1518519
50 6.3481481 0.8195767
51 5.2624339 6.3481481
52 -0.8804233 5.2624339
53 -3.4829976 -0.8804233
54 8.8003358 -3.4829976
55 -2.6993081 8.8003358
56 -2.6326414 -2.6993081
57 0.6840252 -2.6326414
58 -4.7159748 0.6840252
59 -19.6159748 -4.7159748
60 8.1817053 -19.6159748
61 0.5531339 8.1817053
62 -11.3182947 0.5531339
63 2.8959910 -11.3182947
64 -8.4468661 2.8959910
65 -7.1494404 -8.4468661
66 -4.9661070 -7.1494404
67 5.6321123 -4.9661070
68 10.0987790 5.6321123
69 -5.6845543 10.0987790
70 -0.6845543 -5.6845543
71 -12.5845543 -0.6845543
72 -2.6868742 -12.5845543
73 -12.4154457 -2.6868742
74 -3.8868742 -12.4154457
75 12.5274115 -3.8868742
76 9.6845543 12.5274115
77 NA 9.6845543
> dum1 <- dum[2:length(myerror),]
> dum1
lag(myerror, k = 1) myerror
[1,] 3.2853480 -1.3860806
[2,] -16.3860806 3.2853480
[3,] 3.3282051 -16.3860806
[4,] -17.2146520 3.3282051
[5,] -15.2172263 -17.2146520
[6,] -18.6338930 -15.2172263
[7,] -5.1335368 -18.6338930
[8,] 6.8331298 -5.1335368
[9,] -7.9502035 6.8331298
[10,] -9.4502035 -7.9502035
[11,] 9.0497965 -9.4502035
[12,] 9.8474766 9.0497965
[13,] -5.4810948 9.8474766
[14,] 12.4474766 -5.4810948
[15,] -11.4382377 12.4474766
[16,] 4.7189052 -11.4382377
[17,] 7.2163309 4.7189052
[18,] 15.3996642 7.2163309
[19,] -3.9999796 15.3996642
[20,] 0.9666870 -3.9999796
[21,] 10.8833537 0.9666870
[22,] 8.8833537 10.8833537
[23,] 15.4833537 8.8833537
[24,] 12.6810338 15.4833537
[25,] -9.3475376 12.6810338
[26,] 6.3810338 -9.3475376
[27,] -5.3046805 6.3810338
[28,] 8.8524624 -5.3046805
[29,] 9.4498881 8.8524624
[30,] -13.5667786 9.4498881
[31,] -5.5664225 -13.5667786
[32,] -3.5997558 -5.5664225
[33,] -5.6830891 -3.5997558
[34,] 9.5169109 -5.6830891
[35,] 12.4169109 9.5169109
[36,] -15.4854090 12.4169109
[37,] 22.5860195 -15.4854090
[38,] 6.4145910 22.5860195
[39,] -7.2711233 6.4145910
[40,] 3.2860195 -7.2711233
[41,] 9.1834453 3.2860195
[42,] 12.9667786 9.1834453
[43,] 11.7671347 12.9667786
[44,] -11.6661986 11.7671347
[45,] 7.7504681 -11.6661986
[46,] -3.5495319 7.7504681
[47,] -4.7495319 -3.5495319
[48,] -11.1518519 -4.7495319
[49,] 0.8195767 -11.1518519
[50,] 6.3481481 0.8195767
[51,] 5.2624339 6.3481481
[52,] -0.8804233 5.2624339
[53,] -3.4829976 -0.8804233
[54,] 8.8003358 -3.4829976
[55,] -2.6993081 8.8003358
[56,] -2.6326414 -2.6993081
[57,] 0.6840252 -2.6326414
[58,] -4.7159748 0.6840252
[59,] -19.6159748 -4.7159748
[60,] 8.1817053 -19.6159748
[61,] 0.5531339 8.1817053
[62,] -11.3182947 0.5531339
[63,] 2.8959910 -11.3182947
[64,] -8.4468661 2.8959910
[65,] -7.1494404 -8.4468661
[66,] -4.9661070 -7.1494404
[67,] 5.6321123 -4.9661070
[68,] 10.0987790 5.6321123
[69,] -5.6845543 10.0987790
[70,] -0.6845543 -5.6845543
[71,] -12.5845543 -0.6845543
[72,] -2.6868742 -12.5845543
[73,] -12.4154457 -2.6868742
[74,] -3.8868742 -12.4154457
[75,] 12.5274115 -3.8868742
[76,] 9.6845543 12.5274115
> z <- as.data.frame(dum1)
> z
lag(myerror, k = 1) myerror
1 3.2853480 -1.3860806
2 -16.3860806 3.2853480
3 3.3282051 -16.3860806
4 -17.2146520 3.3282051
5 -15.2172263 -17.2146520
6 -18.6338930 -15.2172263
7 -5.1335368 -18.6338930
8 6.8331298 -5.1335368
9 -7.9502035 6.8331298
10 -9.4502035 -7.9502035
11 9.0497965 -9.4502035
12 9.8474766 9.0497965
13 -5.4810948 9.8474766
14 12.4474766 -5.4810948
15 -11.4382377 12.4474766
16 4.7189052 -11.4382377
17 7.2163309 4.7189052
18 15.3996642 7.2163309
19 -3.9999796 15.3996642
20 0.9666870 -3.9999796
21 10.8833537 0.9666870
22 8.8833537 10.8833537
23 15.4833537 8.8833537
24 12.6810338 15.4833537
25 -9.3475376 12.6810338
26 6.3810338 -9.3475376
27 -5.3046805 6.3810338
28 8.8524624 -5.3046805
29 9.4498881 8.8524624
30 -13.5667786 9.4498881
31 -5.5664225 -13.5667786
32 -3.5997558 -5.5664225
33 -5.6830891 -3.5997558
34 9.5169109 -5.6830891
35 12.4169109 9.5169109
36 -15.4854090 12.4169109
37 22.5860195 -15.4854090
38 6.4145910 22.5860195
39 -7.2711233 6.4145910
40 3.2860195 -7.2711233
41 9.1834453 3.2860195
42 12.9667786 9.1834453
43 11.7671347 12.9667786
44 -11.6661986 11.7671347
45 7.7504681 -11.6661986
46 -3.5495319 7.7504681
47 -4.7495319 -3.5495319
48 -11.1518519 -4.7495319
49 0.8195767 -11.1518519
50 6.3481481 0.8195767
51 5.2624339 6.3481481
52 -0.8804233 5.2624339
53 -3.4829976 -0.8804233
54 8.8003358 -3.4829976
55 -2.6993081 8.8003358
56 -2.6326414 -2.6993081
57 0.6840252 -2.6326414
58 -4.7159748 0.6840252
59 -19.6159748 -4.7159748
60 8.1817053 -19.6159748
61 0.5531339 8.1817053
62 -11.3182947 0.5531339
63 2.8959910 -11.3182947
64 -8.4468661 2.8959910
65 -7.1494404 -8.4468661
66 -4.9661070 -7.1494404
67 5.6321123 -4.9661070
68 10.0987790 5.6321123
69 -5.6845543 10.0987790
70 -0.6845543 -5.6845543
71 -12.5845543 -0.6845543
72 -2.6868742 -12.5845543
73 -12.4154457 -2.6868742
74 -3.8868742 -12.4154457
75 12.5274115 -3.8868742
76 9.6845543 12.5274115
> plot(z,main=paste('Residual Lag plot, lowess, and regression line'), ylab='values of Residuals', xlab='lagged values of Residuals')
> lines(lowess(z))
> abline(lm(z))
> grid()
> dev.off()
null device
1
> postscript(file="/var/www/html/rcomp/tmp/77p2i1258731436.ps",horizontal=F,pagecentre=F,paper="special",width=8.3333333333333,height=5.5555555555556)
> acf(mysum$resid, lag.max=length(mysum$resid)/2, main='Residual Autocorrelation Function')
> grid()
> dev.off()
null device
1
> postscript(file="/var/www/html/rcomp/tmp/8x7901258731436.ps",horizontal=F,pagecentre=F,paper="special",width=8.3333333333333,height=5.5555555555556)
> pacf(mysum$resid, lag.max=length(mysum$resid)/2, main='Residual Partial Autocorrelation Function')
> grid()
> dev.off()
null device
1
> postscript(file="/var/www/html/rcomp/tmp/9tr541258731436.ps",horizontal=F,pagecentre=F,paper="special",width=8.3333333333333,height=5.5555555555556)
> opar <- par(mfrow = c(2,2), oma = c(0, 0, 1.1, 0))
> plot(mylm, las = 1, sub='Residual Diagnostics')
> par(opar)
> dev.off()
null device
1
> if (n > n25) {
+ postscript(file="/var/www/html/rcomp/tmp/10bl081258731436.ps",horizontal=F,pagecentre=F,paper="special",width=8.3333333333333,height=5.5555555555556)
+ plot(kp3:nmkm3,gqarr[,2], main='Goldfeld-Quandt test',ylab='2-sided p-value',xlab='breakpoint')
+ grid()
+ 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, 'Multiple Linear Regression - Estimated Regression Equation', 1, TRUE)
> a<-table.row.end(a)
> myeq <- colnames(x)[1]
> myeq <- paste(myeq, '[t] = ', sep='')
> for (i in 1:k){
+ if (mysum$coefficients[i,1] > 0) myeq <- paste(myeq, '+', '')
+ myeq <- paste(myeq, mysum$coefficients[i,1], sep=' ')
+ if (rownames(mysum$coefficients)[i] != '(Intercept)') {
+ myeq <- paste(myeq, rownames(mysum$coefficients)[i], sep='')
+ if (rownames(mysum$coefficients)[i] != 't') myeq <- paste(myeq, '[t]', sep='')
+ }
+ }
> myeq <- paste(myeq, ' + e[t]')
> a<-table.row.start(a)
> a<-table.element(a, myeq)
> a<-table.row.end(a)
> a<-table.end(a)
> table.save(a,file="/var/www/html/rcomp/tmp/11vuvu1258731436.tab")
> a<-table.start()
> a<-table.row.start(a)
> a<-table.element(a,hyperlink('http://www.xycoon.com/ols1.htm','Multiple Linear Regression - Ordinary Least Squares',''), 6, TRUE)
> a<-table.row.end(a)
> a<-table.row.start(a)
> a<-table.element(a,'Variable',header=TRUE)
> a<-table.element(a,'Parameter',header=TRUE)
> a<-table.element(a,'S.D.',header=TRUE)
> a<-table.element(a,'T-STAT
H0: parameter = 0',header=TRUE)
> a<-table.element(a,'2-tail p-value',header=TRUE)
> a<-table.element(a,'1-tail p-value',header=TRUE)
> a<-table.row.end(a)
> for (i in 1:k){
+ a<-table.row.start(a)
+ a<-table.element(a,rownames(mysum$coefficients)[i],header=TRUE)
+ a<-table.element(a,mysum$coefficients[i,1])
+ a<-table.element(a, round(mysum$coefficients[i,2],6))
+ a<-table.element(a, round(mysum$coefficients[i,3],4))
+ a<-table.element(a, round(mysum$coefficients[i,4],6))
+ a<-table.element(a, round(mysum$coefficients[i,4]/2,6))
+ a<-table.row.end(a)
+ }
> a<-table.end(a)
> table.save(a,file="/var/www/html/rcomp/tmp/12dtfk1258731436.tab")
> a<-table.start()
> a<-table.row.start(a)
> a<-table.element(a, 'Multiple Linear Regression - Regression Statistics', 2, TRUE)
> a<-table.row.end(a)
> a<-table.row.start(a)
> a<-table.element(a, 'Multiple R',1,TRUE)
> a<-table.element(a, sqrt(mysum$r.squared))
> a<-table.row.end(a)
> a<-table.row.start(a)
> a<-table.element(a, 'R-squared',1,TRUE)
> a<-table.element(a, mysum$r.squared)
> a<-table.row.end(a)
> a<-table.row.start(a)
> a<-table.element(a, 'Adjusted R-squared',1,TRUE)
> a<-table.element(a, mysum$adj.r.squared)
> a<-table.row.end(a)
> a<-table.row.start(a)
> a<-table.element(a, 'F-TEST (value)',1,TRUE)
> a<-table.element(a, mysum$fstatistic[1])
> a<-table.row.end(a)
> a<-table.row.start(a)
> a<-table.element(a, 'F-TEST (DF numerator)',1,TRUE)
> a<-table.element(a, mysum$fstatistic[2])
> a<-table.row.end(a)
> a<-table.row.start(a)
> a<-table.element(a, 'F-TEST (DF denominator)',1,TRUE)
> a<-table.element(a, mysum$fstatistic[3])
> a<-table.row.end(a)
> a<-table.row.start(a)
> a<-table.element(a, 'p-value',1,TRUE)
> a<-table.element(a, 1-pf(mysum$fstatistic[1],mysum$fstatistic[2],mysum$fstatistic[3]))
> a<-table.row.end(a)
> a<-table.row.start(a)
> a<-table.element(a, 'Multiple Linear Regression - Residual Statistics', 2, TRUE)
> a<-table.row.end(a)
> a<-table.row.start(a)
> a<-table.element(a, 'Residual Standard Deviation',1,TRUE)
> a<-table.element(a, mysum$sigma)
> a<-table.row.end(a)
> a<-table.row.start(a)
> a<-table.element(a, 'Sum Squared Residuals',1,TRUE)
> a<-table.element(a, sum(myerror*myerror))
> a<-table.row.end(a)
> a<-table.end(a)
> table.save(a,file="/var/www/html/rcomp/tmp/13j6ui1258731437.tab")
> a<-table.start()
> a<-table.row.start(a)
> a<-table.element(a, 'Multiple Linear Regression - Actuals, Interpolation, and Residuals', 4, TRUE)
> a<-table.row.end(a)
> a<-table.row.start(a)
> a<-table.element(a, 'Time or Index', 1, TRUE)
> a<-table.element(a, 'Actuals', 1, TRUE)
> a<-table.element(a, 'Interpolation
Forecast', 1, TRUE)
> a<-table.element(a, 'Residuals
Prediction Error', 1, TRUE)
> a<-table.row.end(a)
> for (i in 1:n) {
+ a<-table.row.start(a)
+ a<-table.element(a,i, 1, TRUE)
+ a<-table.element(a,x[i])
+ a<-table.element(a,x[i]-mysum$resid[i])
+ a<-table.element(a,mysum$resid[i])
+ a<-table.row.end(a)
+ }
> a<-table.end(a)
> table.save(a,file="/var/www/html/rcomp/tmp/14oesx1258731437.tab")
> if (n > n25) {
+ a<-table.start()
+ a<-table.row.start(a)
+ a<-table.element(a,'Goldfeld-Quandt test for Heteroskedasticity',4,TRUE)
+ a<-table.row.end(a)
+ a<-table.row.start(a)
+ a<-table.element(a,'p-values',header=TRUE)
+ a<-table.element(a,'Alternative Hypothesis',3,header=TRUE)
+ a<-table.row.end(a)
+ a<-table.row.start(a)
+ a<-table.element(a,'breakpoint index',header=TRUE)
+ a<-table.element(a,'greater',header=TRUE)
+ a<-table.element(a,'2-sided',header=TRUE)
+ a<-table.element(a,'less',header=TRUE)
+ a<-table.row.end(a)
+ for (mypoint in kp3:nmkm3) {
+ a<-table.row.start(a)
+ a<-table.element(a,mypoint,header=TRUE)
+ a<-table.element(a,gqarr[mypoint-kp3+1,1])
+ a<-table.element(a,gqarr[mypoint-kp3+1,2])
+ a<-table.element(a,gqarr[mypoint-kp3+1,3])
+ a<-table.row.end(a)
+ }
+ a<-table.end(a)
+ table.save(a,file="/var/www/html/rcomp/tmp/15fguv1258731437.tab")
+ a<-table.start()
+ a<-table.row.start(a)
+ a<-table.element(a,'Meta Analysis of Goldfeld-Quandt test for Heteroskedasticity',4,TRUE)
+ a<-table.row.end(a)
+ a<-table.row.start(a)
+ a<-table.element(a,'Description',header=TRUE)
+ a<-table.element(a,'# significant tests',header=TRUE)
+ a<-table.element(a,'% significant tests',header=TRUE)
+ a<-table.element(a,'OK/NOK',header=TRUE)
+ a<-table.row.end(a)
+ a<-table.row.start(a)
+ a<-table.element(a,'1% type I error level',header=TRUE)
+ a<-table.element(a,numsignificant1)
+ a<-table.element(a,numsignificant1/numgqtests)
+ if (numsignificant1/numgqtests < 0.01) dum <- 'OK' else dum <- 'NOK'
+ a<-table.element(a,dum)
+ a<-table.row.end(a)
+ a<-table.row.start(a)
+ a<-table.element(a,'5% type I error level',header=TRUE)
+ a<-table.element(a,numsignificant5)
+ a<-table.element(a,numsignificant5/numgqtests)
+ if (numsignificant5/numgqtests < 0.05) dum <- 'OK' else dum <- 'NOK'
+ a<-table.element(a,dum)
+ a<-table.row.end(a)
+ a<-table.row.start(a)
+ a<-table.element(a,'10% type I error level',header=TRUE)
+ a<-table.element(a,numsignificant10)
+ a<-table.element(a,numsignificant10/numgqtests)
+ if (numsignificant10/numgqtests < 0.1) dum <- 'OK' else dum <- 'NOK'
+ a<-table.element(a,dum)
+ a<-table.row.end(a)
+ a<-table.end(a)
+ table.save(a,file="/var/www/html/rcomp/tmp/16bzl41258731437.tab")
+ }
>
> system("convert tmp/1b32o1258731436.ps tmp/1b32o1258731436.png")
> system("convert tmp/2fjcb1258731436.ps tmp/2fjcb1258731436.png")
> system("convert tmp/3hjpd1258731436.ps tmp/3hjpd1258731436.png")
> system("convert tmp/4006s1258731436.ps tmp/4006s1258731436.png")
> system("convert tmp/5d0zx1258731436.ps tmp/5d0zx1258731436.png")
> system("convert tmp/6qvdx1258731436.ps tmp/6qvdx1258731436.png")
> system("convert tmp/77p2i1258731436.ps tmp/77p2i1258731436.png")
> system("convert tmp/8x7901258731436.ps tmp/8x7901258731436.png")
> system("convert tmp/9tr541258731436.ps tmp/9tr541258731436.png")
> system("convert tmp/10bl081258731436.ps tmp/10bl081258731436.png")
>
>
> proc.time()
user system elapsed
2.604 1.573 3.003