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.
> x <- array(list(87.0,106.7,96.3,101.1,107.1,97.8,115.2,113.8,106.1,107.1,89.5,117.5,91.3,113.7,97.6,106.6,100.7,109.8,104.6,108.8,94.7,102.0,101.8,114.5,102.5,116.5,105.3,108.6,110.3,113.9,109.8,109.3,117.3,112.5,118.8,123.4,131.3,115.2,125.9,110.8,133.1,120.4,147.0,117.6,145.8,111.2,164.4,131.1,149.8,118.9,137.7,115.7,151.7,119.6,156.8,113.1,180.0,106.4,180.4,115.5,170.4,111.8,191.6,109.6,199.5,121.5,218.2,109.5,217.5,109.0,205.0,113.4,194.0,112.7,199.3,114.4,219.3,109.2,211.1,116.2,215.2,113.8,240.2,123.6,242.2,112.6,240.7,117.7,255.4,113.3,253.0,110.7,218.2,114.7,203.7,116.9,205.6,120.6,215.6,111.6,188.5,111.9,202.9,116.1,214.0,111.9,230.3,125.1,230.0,115.1,241.0,116.7,259.6,115.8,247.8,116.8,270.3,113.0,289.7,106.5),dim=c(2,60),dimnames=list(c('X','Y'),1:60))
> y <- array(NA,dim=c(2,60),dimnames=list(c('X','Y'),1:60))
> 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
X Y M1 M2 M3 M4 M5 M6 M7 M8 M9 M10 M11 t
1 87.0 106.7 1 0 0 0 0 0 0 0 0 0 0 1
2 96.3 101.1 0 1 0 0 0 0 0 0 0 0 0 2
3 107.1 97.8 0 0 1 0 0 0 0 0 0 0 0 3
4 115.2 113.8 0 0 0 1 0 0 0 0 0 0 0 4
5 106.1 107.1 0 0 0 0 1 0 0 0 0 0 0 5
6 89.5 117.5 0 0 0 0 0 1 0 0 0 0 0 6
7 91.3 113.7 0 0 0 0 0 0 1 0 0 0 0 7
8 97.6 106.6 0 0 0 0 0 0 0 1 0 0 0 8
9 100.7 109.8 0 0 0 0 0 0 0 0 1 0 0 9
10 104.6 108.8 0 0 0 0 0 0 0 0 0 1 0 10
11 94.7 102.0 0 0 0 0 0 0 0 0 0 0 1 11
12 101.8 114.5 0 0 0 0 0 0 0 0 0 0 0 12
13 102.5 116.5 1 0 0 0 0 0 0 0 0 0 0 13
14 105.3 108.6 0 1 0 0 0 0 0 0 0 0 0 14
15 110.3 113.9 0 0 1 0 0 0 0 0 0 0 0 15
16 109.8 109.3 0 0 0 1 0 0 0 0 0 0 0 16
17 117.3 112.5 0 0 0 0 1 0 0 0 0 0 0 17
18 118.8 123.4 0 0 0 0 0 1 0 0 0 0 0 18
19 131.3 115.2 0 0 0 0 0 0 1 0 0 0 0 19
20 125.9 110.8 0 0 0 0 0 0 0 1 0 0 0 20
21 133.1 120.4 0 0 0 0 0 0 0 0 1 0 0 21
22 147.0 117.6 0 0 0 0 0 0 0 0 0 1 0 22
23 145.8 111.2 0 0 0 0 0 0 0 0 0 0 1 23
24 164.4 131.1 0 0 0 0 0 0 0 0 0 0 0 24
25 149.8 118.9 1 0 0 0 0 0 0 0 0 0 0 25
26 137.7 115.7 0 1 0 0 0 0 0 0 0 0 0 26
27 151.7 119.6 0 0 1 0 0 0 0 0 0 0 0 27
28 156.8 113.1 0 0 0 1 0 0 0 0 0 0 0 28
29 180.0 106.4 0 0 0 0 1 0 0 0 0 0 0 29
30 180.4 115.5 0 0 0 0 0 1 0 0 0 0 0 30
31 170.4 111.8 0 0 0 0 0 0 1 0 0 0 0 31
32 191.6 109.6 0 0 0 0 0 0 0 1 0 0 0 32
33 199.5 121.5 0 0 0 0 0 0 0 0 1 0 0 33
34 218.2 109.5 0 0 0 0 0 0 0 0 0 1 0 34
35 217.5 109.0 0 0 0 0 0 0 0 0 0 0 1 35
36 205.0 113.4 0 0 0 0 0 0 0 0 0 0 0 36
37 194.0 112.7 1 0 0 0 0 0 0 0 0 0 0 37
38 199.3 114.4 0 1 0 0 0 0 0 0 0 0 0 38
39 219.3 109.2 0 0 1 0 0 0 0 0 0 0 0 39
40 211.1 116.2 0 0 0 1 0 0 0 0 0 0 0 40
41 215.2 113.8 0 0 0 0 1 0 0 0 0 0 0 41
42 240.2 123.6 0 0 0 0 0 1 0 0 0 0 0 42
43 242.2 112.6 0 0 0 0 0 0 1 0 0 0 0 43
44 240.7 117.7 0 0 0 0 0 0 0 1 0 0 0 44
45 255.4 113.3 0 0 0 0 0 0 0 0 1 0 0 45
46 253.0 110.7 0 0 0 0 0 0 0 0 0 1 0 46
47 218.2 114.7 0 0 0 0 0 0 0 0 0 0 1 47
48 203.7 116.9 0 0 0 0 0 0 0 0 0 0 0 48
49 205.6 120.6 1 0 0 0 0 0 0 0 0 0 0 49
50 215.6 111.6 0 1 0 0 0 0 0 0 0 0 0 50
51 188.5 111.9 0 0 1 0 0 0 0 0 0 0 0 51
52 202.9 116.1 0 0 0 1 0 0 0 0 0 0 0 52
53 214.0 111.9 0 0 0 0 1 0 0 0 0 0 0 53
54 230.3 125.1 0 0 0 0 0 1 0 0 0 0 0 54
55 230.0 115.1 0 0 0 0 0 0 1 0 0 0 0 55
56 241.0 116.7 0 0 0 0 0 0 0 1 0 0 0 56
57 259.6 115.8 0 0 0 0 0 0 0 0 1 0 0 57
58 247.8 116.8 0 0 0 0 0 0 0 0 0 1 0 58
59 270.3 113.0 0 0 0 0 0 0 0 0 0 0 1 59
60 289.7 106.5 0 0 0 0 0 0 0 0 0 0 0 60
> k <- length(x[1,])
> df <- as.data.frame(x)
> (mylm <- lm(df))
Call:
lm(formula = df)
Coefficients:
(Intercept) Y M1 M2 M3 M4
193.402 -0.977 -11.882 -16.659 -15.072 -11.293
M5 M6 M7 M8 M9 M10
-10.364 2.243 -6.876 -5.072 5.871 3.783
M11 t
-6.823 3.148
> (mysum <- summary(mylm))
Call:
lm(formula = df)
Residuals:
Min 1Q Median 3Q Max
-41.041 -12.486 -4.114 12.601 33.107
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 193.4019 65.4495 2.955 0.00492 **
Y -0.9770 0.5732 -1.704 0.09503 .
M1 -11.8816 12.5168 -0.949 0.34745
M2 -16.6592 12.8434 -1.297 0.20106
M3 -15.0716 12.8142 -1.176 0.24558
M4 -11.2934 12.5153 -0.902 0.37156
M5 -10.3641 12.8330 -0.808 0.42347
M6 2.2427 12.7868 0.175 0.86154
M7 -6.8765 12.4989 -0.550 0.58487
M8 -5.0722 12.6104 -0.402 0.68938
M9 5.8709 12.4159 0.473 0.63855
M10 3.7830 12.5811 0.301 0.76501
M11 -6.8228 12.9391 -0.527 0.60052
t 3.1478 0.1606 19.606 < 2e-16 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 19.62 on 46 degrees of freedom
Multiple R-squared: 0.9091, Adjusted R-squared: 0.8834
F-statistic: 35.38 on 13 and 46 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.02662252 0.05324504 0.9733775
[2,] 0.05662729 0.11325459 0.9433727
[3,] 0.10837376 0.21674753 0.8916262
[4,] 0.08199954 0.16399908 0.9180005
[5,] 0.07336908 0.14673817 0.9266309
[6,] 0.08388611 0.16777222 0.9161139
[7,] 0.13691237 0.27382474 0.8630876
[8,] 0.16666502 0.33333004 0.8333350
[9,] 0.14144905 0.28289809 0.8585510
[10,] 0.10596305 0.21192610 0.8940369
[11,] 0.06480090 0.12960179 0.9351991
[12,] 0.04621644 0.09243289 0.9537836
[13,] 0.06358372 0.12716744 0.9364163
[14,] 0.11144619 0.22289237 0.8885538
[15,] 0.12528331 0.25056662 0.8747167
[16,] 0.28569664 0.57139329 0.7143034
[17,] 0.29748596 0.59497191 0.7025140
[18,] 0.35778359 0.71556718 0.6422164
[19,] 0.48253767 0.96507534 0.5174623
[20,] 0.41681499 0.83362998 0.5831850
[21,] 0.62488816 0.75022368 0.3751118
[22,] 0.52679532 0.94640937 0.4732047
[23,] 0.51101357 0.97797285 0.4889864
[24,] 0.46897508 0.93795017 0.5310249
[25,] 0.51850170 0.96299661 0.4814983
[26,] 0.50737260 0.98525479 0.4926274
[27,] 0.45487366 0.90974731 0.5451263
> postscript(file="/var/www/html/rcomp/tmp/1yxpf1227810316.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/2q1r91227810316.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/3ji251227810316.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/4ix8d1227810316.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/52eiq1227810316.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 = 60
Frequency = 1
1 2 3 4 5 6 7
6.580586 12.038965 14.879375 31.685765 11.962499 -10.231069 -6.172397
8 9 10 11 12 13 14
-11.761444 -19.625822 -17.762798 -26.848604 -17.506449 -6.118592 -9.407372
15 16 17 18 19 20 21
-3.964541 -15.884882 -9.335593 -12.940647 -2.480890 -17.131967 -14.643379
22 23 24 25 26 27 28
-4.539002 -4.533997 23.538148 5.752239 -7.844520 5.230475 -2.946215
29 30 31 32 33 34 35
9.630518 3.166818 -4.476809 9.621571 15.057318 20.973058 27.242515
36 37 38 39 40 41 42
9.070761 6.120647 14.711315 24.895376 16.608533 14.286478 33.106695
43 44 45 46 47 48 49
30.330781 28.861449 25.171675 19.171458 -4.262469 -26.583680 -12.334880
50 51 52 53 54 55 56
-9.498388 -41.040686 -29.463200 -26.543902 -13.101797 -17.200686 -9.589608
57 58 59 60
-5.959792 -17.842716 8.402556 11.481220
> postscript(file="/var/www/html/rcomp/tmp/6lcg71227810316.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 = 60
Frequency = 1
lag(myerror, k = 1) myerror
0 6.580586 NA
1 12.038965 6.580586
2 14.879375 12.038965
3 31.685765 14.879375
4 11.962499 31.685765
5 -10.231069 11.962499
6 -6.172397 -10.231069
7 -11.761444 -6.172397
8 -19.625822 -11.761444
9 -17.762798 -19.625822
10 -26.848604 -17.762798
11 -17.506449 -26.848604
12 -6.118592 -17.506449
13 -9.407372 -6.118592
14 -3.964541 -9.407372
15 -15.884882 -3.964541
16 -9.335593 -15.884882
17 -12.940647 -9.335593
18 -2.480890 -12.940647
19 -17.131967 -2.480890
20 -14.643379 -17.131967
21 -4.539002 -14.643379
22 -4.533997 -4.539002
23 23.538148 -4.533997
24 5.752239 23.538148
25 -7.844520 5.752239
26 5.230475 -7.844520
27 -2.946215 5.230475
28 9.630518 -2.946215
29 3.166818 9.630518
30 -4.476809 3.166818
31 9.621571 -4.476809
32 15.057318 9.621571
33 20.973058 15.057318
34 27.242515 20.973058
35 9.070761 27.242515
36 6.120647 9.070761
37 14.711315 6.120647
38 24.895376 14.711315
39 16.608533 24.895376
40 14.286478 16.608533
41 33.106695 14.286478
42 30.330781 33.106695
43 28.861449 30.330781
44 25.171675 28.861449
45 19.171458 25.171675
46 -4.262469 19.171458
47 -26.583680 -4.262469
48 -12.334880 -26.583680
49 -9.498388 -12.334880
50 -41.040686 -9.498388
51 -29.463200 -41.040686
52 -26.543902 -29.463200
53 -13.101797 -26.543902
54 -17.200686 -13.101797
55 -9.589608 -17.200686
56 -5.959792 -9.589608
57 -17.842716 -5.959792
58 8.402556 -17.842716
59 11.481220 8.402556
60 NA 11.481220
> dum1 <- dum[2:length(myerror),]
> dum1
lag(myerror, k = 1) myerror
[1,] 12.038965 6.580586
[2,] 14.879375 12.038965
[3,] 31.685765 14.879375
[4,] 11.962499 31.685765
[5,] -10.231069 11.962499
[6,] -6.172397 -10.231069
[7,] -11.761444 -6.172397
[8,] -19.625822 -11.761444
[9,] -17.762798 -19.625822
[10,] -26.848604 -17.762798
[11,] -17.506449 -26.848604
[12,] -6.118592 -17.506449
[13,] -9.407372 -6.118592
[14,] -3.964541 -9.407372
[15,] -15.884882 -3.964541
[16,] -9.335593 -15.884882
[17,] -12.940647 -9.335593
[18,] -2.480890 -12.940647
[19,] -17.131967 -2.480890
[20,] -14.643379 -17.131967
[21,] -4.539002 -14.643379
[22,] -4.533997 -4.539002
[23,] 23.538148 -4.533997
[24,] 5.752239 23.538148
[25,] -7.844520 5.752239
[26,] 5.230475 -7.844520
[27,] -2.946215 5.230475
[28,] 9.630518 -2.946215
[29,] 3.166818 9.630518
[30,] -4.476809 3.166818
[31,] 9.621571 -4.476809
[32,] 15.057318 9.621571
[33,] 20.973058 15.057318
[34,] 27.242515 20.973058
[35,] 9.070761 27.242515
[36,] 6.120647 9.070761
[37,] 14.711315 6.120647
[38,] 24.895376 14.711315
[39,] 16.608533 24.895376
[40,] 14.286478 16.608533
[41,] 33.106695 14.286478
[42,] 30.330781 33.106695
[43,] 28.861449 30.330781
[44,] 25.171675 28.861449
[45,] 19.171458 25.171675
[46,] -4.262469 19.171458
[47,] -26.583680 -4.262469
[48,] -12.334880 -26.583680
[49,] -9.498388 -12.334880
[50,] -41.040686 -9.498388
[51,] -29.463200 -41.040686
[52,] -26.543902 -29.463200
[53,] -13.101797 -26.543902
[54,] -17.200686 -13.101797
[55,] -9.589608 -17.200686
[56,] -5.959792 -9.589608
[57,] -17.842716 -5.959792
[58,] 8.402556 -17.842716
[59,] 11.481220 8.402556
> z <- as.data.frame(dum1)
> z
lag(myerror, k = 1) myerror
1 12.038965 6.580586
2 14.879375 12.038965
3 31.685765 14.879375
4 11.962499 31.685765
5 -10.231069 11.962499
6 -6.172397 -10.231069
7 -11.761444 -6.172397
8 -19.625822 -11.761444
9 -17.762798 -19.625822
10 -26.848604 -17.762798
11 -17.506449 -26.848604
12 -6.118592 -17.506449
13 -9.407372 -6.118592
14 -3.964541 -9.407372
15 -15.884882 -3.964541
16 -9.335593 -15.884882
17 -12.940647 -9.335593
18 -2.480890 -12.940647
19 -17.131967 -2.480890
20 -14.643379 -17.131967
21 -4.539002 -14.643379
22 -4.533997 -4.539002
23 23.538148 -4.533997
24 5.752239 23.538148
25 -7.844520 5.752239
26 5.230475 -7.844520
27 -2.946215 5.230475
28 9.630518 -2.946215
29 3.166818 9.630518
30 -4.476809 3.166818
31 9.621571 -4.476809
32 15.057318 9.621571
33 20.973058 15.057318
34 27.242515 20.973058
35 9.070761 27.242515
36 6.120647 9.070761
37 14.711315 6.120647
38 24.895376 14.711315
39 16.608533 24.895376
40 14.286478 16.608533
41 33.106695 14.286478
42 30.330781 33.106695
43 28.861449 30.330781
44 25.171675 28.861449
45 19.171458 25.171675
46 -4.262469 19.171458
47 -26.583680 -4.262469
48 -12.334880 -26.583680
49 -9.498388 -12.334880
50 -41.040686 -9.498388
51 -29.463200 -41.040686
52 -26.543902 -29.463200
53 -13.101797 -26.543902
54 -17.200686 -13.101797
55 -9.589608 -17.200686
56 -5.959792 -9.589608
57 -17.842716 -5.959792
58 8.402556 -17.842716
59 11.481220 8.402556
> 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/76mdc1227810316.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/8hxlk1227810316.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/9mt8r1227810316.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/10wt9p1227810316.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/116bab1227810317.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/12g5au1227810317.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/1303vb1227810317.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/14j4x71227810317.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/15fpje1227810317.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/16xvof1227810317.tab")
+ }
>
> system("convert tmp/1yxpf1227810316.ps tmp/1yxpf1227810316.png")
> system("convert tmp/2q1r91227810316.ps tmp/2q1r91227810316.png")
> system("convert tmp/3ji251227810316.ps tmp/3ji251227810316.png")
> system("convert tmp/4ix8d1227810316.ps tmp/4ix8d1227810316.png")
> system("convert tmp/52eiq1227810316.ps tmp/52eiq1227810316.png")
> system("convert tmp/6lcg71227810316.ps tmp/6lcg71227810316.png")
> system("convert tmp/76mdc1227810316.ps tmp/76mdc1227810316.png")
> system("convert tmp/8hxlk1227810316.ps tmp/8hxlk1227810316.png")
> system("convert tmp/9mt8r1227810316.ps tmp/9mt8r1227810316.png")
> system("convert tmp/10wt9p1227810316.ps tmp/10wt9p1227810316.png")
>
>
> proc.time()
user system elapsed
2.388 1.543 2.778