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(-22,46,-20,50,-17,49,-21,48,-16,50,-11,47,-19,50,-31,49,-36,51,-33,52,-26,48,-38,55,-27,56,-21,43,-17,44,-14,50,-16,49,-16,47,-15,46,-7,50,-9,49,2,53,-6,54,0,56,7,56,4,58,-5,53,2,51,0,52,3,53,10,56,4,54,5,54,7,56,1,59,-8,62,-3,62,-16,73,-22,76,-32,80,-30,77,-32,81,-38,80,-41,80,-46,81,-58,80,-55,77,-48,71,-58,71,-58,64,-68,64,-75,47,-77,41,-75,35,-71,34,-63,33,-61,23,-53,16,-41,16,-35,8,-33,9),dim=c(2,61),dimnames=list(c('Econo','Price'),1:61))
> y <- array(NA,dim=c(2,61),dimnames=list(c('Econo','Price'),1:61))
> for (i in 1:dim(x)[1])
+ {
+ for (j in 1:dim(x)[2])
+ {
+ y[i,j] <- as.numeric(x[i,j])
+ }
+ }
> par3 = 'No 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
Econo Price M1 M2 M3 M4 M5 M6 M7 M8 M9 M10 M11
1 -22 46 1 0 0 0 0 0 0 0 0 0 0
2 -20 50 0 1 0 0 0 0 0 0 0 0 0
3 -17 49 0 0 1 0 0 0 0 0 0 0 0
4 -21 48 0 0 0 1 0 0 0 0 0 0 0
5 -16 50 0 0 0 0 1 0 0 0 0 0 0
6 -11 47 0 0 0 0 0 1 0 0 0 0 0
7 -19 50 0 0 0 0 0 0 1 0 0 0 0
8 -31 49 0 0 0 0 0 0 0 1 0 0 0
9 -36 51 0 0 0 0 0 0 0 0 1 0 0
10 -33 52 0 0 0 0 0 0 0 0 0 1 0
11 -26 48 0 0 0 0 0 0 0 0 0 0 1
12 -38 55 0 0 0 0 0 0 0 0 0 0 0
13 -27 56 1 0 0 0 0 0 0 0 0 0 0
14 -21 43 0 1 0 0 0 0 0 0 0 0 0
15 -17 44 0 0 1 0 0 0 0 0 0 0 0
16 -14 50 0 0 0 1 0 0 0 0 0 0 0
17 -16 49 0 0 0 0 1 0 0 0 0 0 0
18 -16 47 0 0 0 0 0 1 0 0 0 0 0
19 -15 46 0 0 0 0 0 0 1 0 0 0 0
20 -7 50 0 0 0 0 0 0 0 1 0 0 0
21 -9 49 0 0 0 0 0 0 0 0 1 0 0
22 2 53 0 0 0 0 0 0 0 0 0 1 0
23 -6 54 0 0 0 0 0 0 0 0 0 0 1
24 0 56 0 0 0 0 0 0 0 0 0 0 0
25 7 56 1 0 0 0 0 0 0 0 0 0 0
26 4 58 0 1 0 0 0 0 0 0 0 0 0
27 -5 53 0 0 1 0 0 0 0 0 0 0 0
28 2 51 0 0 0 1 0 0 0 0 0 0 0
29 0 52 0 0 0 0 1 0 0 0 0 0 0
30 3 53 0 0 0 0 0 1 0 0 0 0 0
31 10 56 0 0 0 0 0 0 1 0 0 0 0
32 4 54 0 0 0 0 0 0 0 1 0 0 0
33 5 54 0 0 0 0 0 0 0 0 1 0 0
34 7 56 0 0 0 0 0 0 0 0 0 1 0
35 1 59 0 0 0 0 0 0 0 0 0 0 1
36 -8 62 0 0 0 0 0 0 0 0 0 0 0
37 -3 62 1 0 0 0 0 0 0 0 0 0 0
38 -16 73 0 1 0 0 0 0 0 0 0 0 0
39 -22 76 0 0 1 0 0 0 0 0 0 0 0
40 -32 80 0 0 0 1 0 0 0 0 0 0 0
41 -30 77 0 0 0 0 1 0 0 0 0 0 0
42 -32 81 0 0 0 0 0 1 0 0 0 0 0
43 -38 80 0 0 0 0 0 0 1 0 0 0 0
44 -41 80 0 0 0 0 0 0 0 1 0 0 0
45 -46 81 0 0 0 0 0 0 0 0 1 0 0
46 -58 80 0 0 0 0 0 0 0 0 0 1 0
47 -55 77 0 0 0 0 0 0 0 0 0 0 1
48 -48 71 0 0 0 0 0 0 0 0 0 0 0
49 -58 71 1 0 0 0 0 0 0 0 0 0 0
50 -58 64 0 1 0 0 0 0 0 0 0 0 0
51 -68 64 0 0 1 0 0 0 0 0 0 0 0
52 -75 47 0 0 0 1 0 0 0 0 0 0 0
53 -77 41 0 0 0 0 1 0 0 0 0 0 0
54 -75 35 0 0 0 0 0 1 0 0 0 0 0
55 -71 34 0 0 0 0 0 0 1 0 0 0 0
56 -63 33 0 0 0 0 0 0 0 1 0 0 0
57 -61 23 0 0 0 0 0 0 0 0 1 0 0
58 -53 16 0 0 0 0 0 0 0 0 0 1 0
59 -41 16 0 0 0 0 0 0 0 0 0 0 1
60 -35 8 0 0 0 0 0 0 0 0 0 0 0
61 -33 9 1 0 0 0 0 0 0 0 0 0 0
> k <- length(x[1,])
> df <- as.data.frame(x)
> (mylm <- lm(df))
Call:
lm(formula = df)
Coefficients:
(Intercept) Price M1 M2 M3 M4
-32.6420 0.1358 3.1876 2.6226 -0.9231 -2.8516
M5 M6 M7 M8 M9 M10
-2.4616 -0.6987 -1.1801 -2.1801 -3.7629 -1.3358
M11
0.3457
> (mysum <- summary(mylm))
Call:
lm(formula = df)
Residuals:
Min 1Q Median 3Q Max
-47.462 -15.038 3.182 18.966 36.220
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) -32.6420 15.8199 -2.063 0.0445 *
Price 0.1358 0.2061 0.659 0.5133
M1 3.1876 16.1557 0.197 0.8444
M2 2.6226 16.9390 0.155 0.8776
M3 -0.9231 16.9320 -0.055 0.9567
M4 -2.8516 16.9029 -0.169 0.8667
M5 -2.4616 16.8884 -0.146 0.8847
M6 -0.6987 16.8800 -0.041 0.9672
M7 -1.1801 16.8837 -0.070 0.9446
M8 -2.1801 16.8837 -0.129 0.8978
M9 -3.7629 16.8757 -0.223 0.8245
M10 -1.3358 16.8751 -0.079 0.9372
M11 0.3457 16.8741 0.020 0.9837
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 26.68 on 48 degrees of freedom
Multiple R-squared: 0.01619, Adjusted R-squared: -0.2298
F-statistic: 0.06584 on 12 and 48 DF, p-value: 1
> 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,] 1.745089e-03 3.490178e-03 0.99825491
[2,] 1.524120e-04 3.048240e-04 0.99984759
[3,] 2.546852e-05 5.093704e-05 0.99997453
[4,] 2.886651e-06 5.773302e-06 0.99999711
[5,] 1.512160e-04 3.024321e-04 0.99984878
[6,] 4.249668e-04 8.499337e-04 0.99957503
[7,] 1.590599e-03 3.181197e-03 0.99840940
[8,] 9.783756e-04 1.956751e-03 0.99902162
[9,] 2.048227e-03 4.096455e-03 0.99795177
[10,] 2.359832e-03 4.719665e-03 0.99764017
[11,] 1.275581e-03 2.551163e-03 0.99872442
[12,] 6.462014e-04 1.292403e-03 0.99935380
[13,] 6.085653e-04 1.217131e-03 0.99939143
[14,] 4.995811e-04 9.991621e-04 0.99950042
[15,] 4.191537e-04 8.383074e-04 0.99958085
[16,] 5.797888e-04 1.159578e-03 0.99942021
[17,] 8.089584e-04 1.617917e-03 0.99919104
[18,] 1.757720e-03 3.515439e-03 0.99824228
[19,] 4.631867e-03 9.263735e-03 0.99536813
[20,] 6.319078e-03 1.263816e-02 0.99368092
[21,] 4.855995e-03 9.711991e-03 0.99514400
[22,] 6.485095e-03 1.297019e-02 0.99351490
[23,] 3.131242e-02 6.262483e-02 0.96868758
[24,] 8.035363e-02 1.607073e-01 0.91964637
[25,] 1.470257e-01 2.940514e-01 0.85297432
[26,] 2.672108e-01 5.344215e-01 0.73278924
[27,] 4.672139e-01 9.344278e-01 0.53278611
[28,] 6.683399e-01 6.633201e-01 0.33166006
[29,] 7.833482e-01 4.333036e-01 0.21665181
[30,] 9.345247e-01 1.309506e-01 0.06547528
> postscript(file="/var/www/html/rcomp/tmp/19fg81258815334.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/266y71258815334.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/3f83o1258815334.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/4u4vv1258815334.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/5fi941258815334.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 = 61
Frequency = 1
1 2 3 4 5 6
1.2096820 3.2317292 9.9131815 7.9774277 12.3158646 15.9602215
7 8 9 10 11 12
8.0344123 -2.8298339 -6.5185477 -6.0814523 -0.2198892 -12.8244677
13 14 15 16 17 18
-5.1478564 3.1820061 10.5919507 14.7059200 12.4516184 10.9602215
19 20 21 22 23 24
12.5774277 21.0344123 20.7529600 28.7827939 18.9655877 25.0397785
25 26 27 28 29 30
28.8521436 26.1456985 21.3701661 30.5701661 28.0443569 29.1456985
31 32 33 34 35 36
36.2198892 31.4913969 34.0741908 33.3755323 25.2868185 16.2252554
37 38 39 40 41 42
18.0376205 4.1093908 1.2478277 -7.3666953 -5.3494892 -9.6554091
43 44 45 46 47 48
-15.0382030 -17.0382030 -20.5911630 -34.8825599 -33.1567507 -24.9965292
49 50 51 52 53 54
-38.1841640 -36.6688246 -43.1231261 -45.8868185 -47.4623508 -46.4107324
55 56 57 58 59 60
-41.7935262 -32.6577724 -27.7174401 -21.1943140 -10.8757663 -3.4440370
61
-4.7674258
> postscript(file="/var/www/html/rcomp/tmp/6i8tm1258815334.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 = 61
Frequency = 1
lag(myerror, k = 1) myerror
0 1.2096820 NA
1 3.2317292 1.2096820
2 9.9131815 3.2317292
3 7.9774277 9.9131815
4 12.3158646 7.9774277
5 15.9602215 12.3158646
6 8.0344123 15.9602215
7 -2.8298339 8.0344123
8 -6.5185477 -2.8298339
9 -6.0814523 -6.5185477
10 -0.2198892 -6.0814523
11 -12.8244677 -0.2198892
12 -5.1478564 -12.8244677
13 3.1820061 -5.1478564
14 10.5919507 3.1820061
15 14.7059200 10.5919507
16 12.4516184 14.7059200
17 10.9602215 12.4516184
18 12.5774277 10.9602215
19 21.0344123 12.5774277
20 20.7529600 21.0344123
21 28.7827939 20.7529600
22 18.9655877 28.7827939
23 25.0397785 18.9655877
24 28.8521436 25.0397785
25 26.1456985 28.8521436
26 21.3701661 26.1456985
27 30.5701661 21.3701661
28 28.0443569 30.5701661
29 29.1456985 28.0443569
30 36.2198892 29.1456985
31 31.4913969 36.2198892
32 34.0741908 31.4913969
33 33.3755323 34.0741908
34 25.2868185 33.3755323
35 16.2252554 25.2868185
36 18.0376205 16.2252554
37 4.1093908 18.0376205
38 1.2478277 4.1093908
39 -7.3666953 1.2478277
40 -5.3494892 -7.3666953
41 -9.6554091 -5.3494892
42 -15.0382030 -9.6554091
43 -17.0382030 -15.0382030
44 -20.5911630 -17.0382030
45 -34.8825599 -20.5911630
46 -33.1567507 -34.8825599
47 -24.9965292 -33.1567507
48 -38.1841640 -24.9965292
49 -36.6688246 -38.1841640
50 -43.1231261 -36.6688246
51 -45.8868185 -43.1231261
52 -47.4623508 -45.8868185
53 -46.4107324 -47.4623508
54 -41.7935262 -46.4107324
55 -32.6577724 -41.7935262
56 -27.7174401 -32.6577724
57 -21.1943140 -27.7174401
58 -10.8757663 -21.1943140
59 -3.4440370 -10.8757663
60 -4.7674258 -3.4440370
61 NA -4.7674258
> dum1 <- dum[2:length(myerror),]
> dum1
lag(myerror, k = 1) myerror
[1,] 3.2317292 1.2096820
[2,] 9.9131815 3.2317292
[3,] 7.9774277 9.9131815
[4,] 12.3158646 7.9774277
[5,] 15.9602215 12.3158646
[6,] 8.0344123 15.9602215
[7,] -2.8298339 8.0344123
[8,] -6.5185477 -2.8298339
[9,] -6.0814523 -6.5185477
[10,] -0.2198892 -6.0814523
[11,] -12.8244677 -0.2198892
[12,] -5.1478564 -12.8244677
[13,] 3.1820061 -5.1478564
[14,] 10.5919507 3.1820061
[15,] 14.7059200 10.5919507
[16,] 12.4516184 14.7059200
[17,] 10.9602215 12.4516184
[18,] 12.5774277 10.9602215
[19,] 21.0344123 12.5774277
[20,] 20.7529600 21.0344123
[21,] 28.7827939 20.7529600
[22,] 18.9655877 28.7827939
[23,] 25.0397785 18.9655877
[24,] 28.8521436 25.0397785
[25,] 26.1456985 28.8521436
[26,] 21.3701661 26.1456985
[27,] 30.5701661 21.3701661
[28,] 28.0443569 30.5701661
[29,] 29.1456985 28.0443569
[30,] 36.2198892 29.1456985
[31,] 31.4913969 36.2198892
[32,] 34.0741908 31.4913969
[33,] 33.3755323 34.0741908
[34,] 25.2868185 33.3755323
[35,] 16.2252554 25.2868185
[36,] 18.0376205 16.2252554
[37,] 4.1093908 18.0376205
[38,] 1.2478277 4.1093908
[39,] -7.3666953 1.2478277
[40,] -5.3494892 -7.3666953
[41,] -9.6554091 -5.3494892
[42,] -15.0382030 -9.6554091
[43,] -17.0382030 -15.0382030
[44,] -20.5911630 -17.0382030
[45,] -34.8825599 -20.5911630
[46,] -33.1567507 -34.8825599
[47,] -24.9965292 -33.1567507
[48,] -38.1841640 -24.9965292
[49,] -36.6688246 -38.1841640
[50,] -43.1231261 -36.6688246
[51,] -45.8868185 -43.1231261
[52,] -47.4623508 -45.8868185
[53,] -46.4107324 -47.4623508
[54,] -41.7935262 -46.4107324
[55,] -32.6577724 -41.7935262
[56,] -27.7174401 -32.6577724
[57,] -21.1943140 -27.7174401
[58,] -10.8757663 -21.1943140
[59,] -3.4440370 -10.8757663
[60,] -4.7674258 -3.4440370
> z <- as.data.frame(dum1)
> z
lag(myerror, k = 1) myerror
1 3.2317292 1.2096820
2 9.9131815 3.2317292
3 7.9774277 9.9131815
4 12.3158646 7.9774277
5 15.9602215 12.3158646
6 8.0344123 15.9602215
7 -2.8298339 8.0344123
8 -6.5185477 -2.8298339
9 -6.0814523 -6.5185477
10 -0.2198892 -6.0814523
11 -12.8244677 -0.2198892
12 -5.1478564 -12.8244677
13 3.1820061 -5.1478564
14 10.5919507 3.1820061
15 14.7059200 10.5919507
16 12.4516184 14.7059200
17 10.9602215 12.4516184
18 12.5774277 10.9602215
19 21.0344123 12.5774277
20 20.7529600 21.0344123
21 28.7827939 20.7529600
22 18.9655877 28.7827939
23 25.0397785 18.9655877
24 28.8521436 25.0397785
25 26.1456985 28.8521436
26 21.3701661 26.1456985
27 30.5701661 21.3701661
28 28.0443569 30.5701661
29 29.1456985 28.0443569
30 36.2198892 29.1456985
31 31.4913969 36.2198892
32 34.0741908 31.4913969
33 33.3755323 34.0741908
34 25.2868185 33.3755323
35 16.2252554 25.2868185
36 18.0376205 16.2252554
37 4.1093908 18.0376205
38 1.2478277 4.1093908
39 -7.3666953 1.2478277
40 -5.3494892 -7.3666953
41 -9.6554091 -5.3494892
42 -15.0382030 -9.6554091
43 -17.0382030 -15.0382030
44 -20.5911630 -17.0382030
45 -34.8825599 -20.5911630
46 -33.1567507 -34.8825599
47 -24.9965292 -33.1567507
48 -38.1841640 -24.9965292
49 -36.6688246 -38.1841640
50 -43.1231261 -36.6688246
51 -45.8868185 -43.1231261
52 -47.4623508 -45.8868185
53 -46.4107324 -47.4623508
54 -41.7935262 -46.4107324
55 -32.6577724 -41.7935262
56 -27.7174401 -32.6577724
57 -21.1943140 -27.7174401
58 -10.8757663 -21.1943140
59 -3.4440370 -10.8757663
60 -4.7674258 -3.4440370
> 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/75zlz1258815334.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/8ikzn1258815334.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/9tlvv1258815334.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/10ek441258815334.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/114a6m1258815335.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/12jndx1258815335.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/131lr81258815335.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/14p2ab1258815335.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/15gq6v1258815335.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/165fkh1258815335.tab")
+ }
>
> system("convert tmp/19fg81258815334.ps tmp/19fg81258815334.png")
> system("convert tmp/266y71258815334.ps tmp/266y71258815334.png")
> system("convert tmp/3f83o1258815334.ps tmp/3f83o1258815334.png")
> system("convert tmp/4u4vv1258815334.ps tmp/4u4vv1258815334.png")
> system("convert tmp/5fi941258815334.ps tmp/5fi941258815334.png")
> system("convert tmp/6i8tm1258815334.ps tmp/6i8tm1258815334.png")
> system("convert tmp/75zlz1258815334.ps tmp/75zlz1258815334.png")
> system("convert tmp/8ikzn1258815334.ps tmp/8ikzn1258815334.png")
> system("convert tmp/9tlvv1258815334.ps tmp/9tlvv1258815334.png")
> system("convert tmp/10ek441258815334.ps tmp/10ek441258815334.png")
>
>
> proc.time()
user system elapsed
2.412 1.563 3.320