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(3484.74
+ ,13830.14
+ ,9349.44
+ ,7977
+ ,-5.6
+ ,6
+ ,1
+ ,2.77
+ ,3411.13
+ ,14153.22
+ ,9327.78
+ ,8241
+ ,-6.2
+ ,3
+ ,1
+ ,2.76
+ ,3288.18
+ ,15418.03
+ ,9753.63
+ ,8444
+ ,-7.1
+ ,2
+ ,1.2
+ ,2.76
+ ,3280.37
+ ,16666.97
+ ,10443.5
+ ,8490
+ ,-1.4
+ ,2
+ ,1.2
+ ,2.46
+ ,3173.95
+ ,16505.21
+ ,10853.87
+ ,8388
+ ,-0.1
+ ,2
+ ,0.8
+ ,2.46
+ ,3165.26
+ ,17135.96
+ ,10704.02
+ ,8099
+ ,-0.9
+ ,-8
+ ,0.7
+ ,2.47
+ ,3092.71
+ ,18033.25
+ ,11052.23
+ ,7984
+ ,0
+ ,0
+ ,0.7
+ ,2.71
+ ,3053.05
+ ,17671
+ ,10935.47
+ ,7786
+ ,0.1
+ ,-2
+ ,0.9
+ ,2.8
+ ,3181.96
+ ,17544.22
+ ,10714.03
+ ,8086
+ ,2.6
+ ,3
+ ,1.2
+ ,2.89
+ ,2999.93
+ ,17677.9
+ ,10394.48
+ ,9315
+ ,6
+ ,5
+ ,1.3
+ ,3.36
+ ,3249.57
+ ,18470.97
+ ,10817.9
+ ,9113
+ ,6.4
+ ,8
+ ,1.5
+ ,3.31
+ ,3210.52
+ ,18409.96
+ ,11251.2
+ ,9023
+ ,8.6
+ ,8
+ ,1.9
+ ,3.5
+ ,3030.29
+ ,18941.6
+ ,11281.26
+ ,9026
+ ,6.4
+ ,9
+ ,1.8
+ ,3.51
+ ,2803.47
+ ,19685.53
+ ,10539.68
+ ,9787
+ ,7.7
+ ,11
+ ,1.9
+ ,3.71
+ ,2767.63
+ ,19834.71
+ ,10483.39
+ ,9536
+ ,9.2
+ ,13
+ ,2.2
+ ,3.71
+ ,2882.6
+ ,19598.93
+ ,10947.43
+ ,9490
+ ,8.6
+ ,12
+ ,2.1
+ ,3.71
+ ,2863.36
+ ,17039.97
+ ,10580.27
+ ,9736
+ ,7.4
+ ,13
+ ,2.2
+ ,4.21
+ ,2897.06
+ ,16969.28
+ ,10582.92
+ ,9694
+ ,8.6
+ ,15
+ ,2.7
+ ,4.21
+ ,3012.61
+ ,16973.38
+ ,10654.41
+ ,9647
+ ,6.2
+ ,13
+ ,2.8
+ ,4.21
+ ,3142.95
+ ,16329.89
+ ,11014.51
+ ,9753
+ ,6
+ ,16
+ ,2.9
+ ,4.5
+ ,3032.93
+ ,16153.34
+ ,10967.87
+ ,10070
+ ,6.6
+ ,10
+ ,3.4
+ ,4.51
+ ,3045.78
+ ,15311.7
+ ,10433.56
+ ,10137
+ ,5.1
+ ,14
+ ,3
+ ,4.51
+ ,3110.52
+ ,14760.87
+ ,10665.78
+ ,9984
+ ,4.7
+ ,14
+ ,3.1
+ ,4.51
+ ,3013.24
+ ,14452.93
+ ,10666.71
+ ,9732
+ ,5
+ ,15
+ ,2.5
+ ,4.32
+ ,2987.1
+ ,13720.95
+ ,10682.74
+ ,9103
+ ,3.6
+ ,13
+ ,2.2
+ ,4.02
+ ,2995.55
+ ,13266.27
+ ,10777.22
+ ,9155
+ ,1.9
+ ,8
+ ,2.3
+ ,4.02
+ ,2833.18
+ ,12708.47
+ ,10052.6
+ ,9308
+ ,-0.1
+ ,7
+ ,2.1
+ ,3.85
+ ,2848.96
+ ,13411.84
+ ,10213.97
+ ,9394
+ ,-5.7
+ ,3
+ ,2.8
+ ,3.84
+ ,2794.83
+ ,13975.55
+ ,10546.82
+ ,9948
+ ,-5.6
+ ,3
+ ,3.1
+ ,4.02
+ ,2845.26
+ ,12974.89
+ ,10767.2
+ ,10177
+ ,-6.4
+ ,4
+ ,2.9
+ ,3.82
+ ,2915.02
+ ,12151.11
+ ,10444.5
+ ,10002
+ ,-7.7
+ ,4
+ ,2.6
+ ,3.75
+ ,2892.63
+ ,11576.21
+ ,10314.68
+ ,9728
+ ,-8
+ ,0
+ ,2.7
+ ,3.74
+ ,2604.42
+ ,9996.83
+ ,9042.56
+ ,10002
+ ,-11.9
+ ,-4
+ ,2.3
+ ,3.14
+ ,2641.65
+ ,10438.9
+ ,9220.75
+ ,10063
+ ,-15.4
+ ,-14
+ ,2.3
+ ,2.91
+ ,2659.81
+ ,10511.22
+ ,9721.84
+ ,10018
+ ,-15.5
+ ,-18
+ ,2.1
+ ,2.84
+ ,2638.53
+ ,10496.2
+ ,9978.53
+ ,9960
+ ,-13.4
+ ,-8
+ ,2.2
+ ,2.85
+ ,2720.25
+ ,10300.79
+ ,9923.81
+ ,10236
+ ,-10.9
+ ,-1
+ ,2.9
+ ,2.85
+ ,2745.88
+ ,9981.65
+ ,9892.56
+ ,10893
+ ,-10.8
+ ,1
+ ,2.6
+ ,3.08
+ ,2735.7
+ ,11448.79
+ ,10500.98
+ ,10756
+ ,-7.3
+ ,2
+ ,2.7
+ ,3.3
+ ,2811.7
+ ,11384.49
+ ,10179.35
+ ,10940
+ ,-6.5
+ ,0
+ ,1.8
+ ,3.29
+ ,2799.43
+ ,11717.46
+ ,10080.48
+ ,10997
+ ,-5.1
+ ,1
+ ,1.3
+ ,3.26
+ ,2555.28
+ ,10965.88
+ ,9492.44
+ ,10827
+ ,-5.3
+ ,0
+ ,0.9
+ ,3.26
+ ,2304.98
+ ,10352.27
+ ,8616.49
+ ,10166
+ ,-6.8
+ ,-1
+ ,1.3
+ ,3.11
+ ,2214.95
+ ,9751.2
+ ,8685.4
+ ,10186
+ ,-8.4
+ ,-3
+ ,1.3
+ ,2.84
+ ,2065.81
+ ,9354.01
+ ,8160.67
+ ,10457
+ ,-8.4
+ ,-3
+ ,1.3
+ ,2.71
+ ,1940.49
+ ,8792.5
+ ,8048.1
+ ,10368
+ ,-9.7
+ ,-3
+ ,1.3
+ ,2.69
+ ,2042
+ ,8721.14
+ ,8641.21
+ ,10244
+ ,-8.8
+ ,-4
+ ,1.1
+ ,2.65
+ ,1995.37
+ ,8692.94
+ ,8526.63
+ ,10511
+ ,-9.6
+ ,-8
+ ,1.4
+ ,2.57
+ ,1946.81
+ ,8570.73
+ ,8474.21
+ ,10812
+ ,-11.5
+ ,-9
+ ,1.2
+ ,2.32
+ ,1765.9
+ ,8538.47
+ ,7916.13
+ ,10738
+ ,-11
+ ,-13
+ ,1.7
+ ,2.12
+ ,1635.25
+ ,8169.75
+ ,7977.64
+ ,10171
+ ,-14.9
+ ,-18
+ ,1.8
+ ,2.05)
+ ,dim=c(8
+ ,51)
+ ,dimnames=list(c('BEL_20'
+ ,'Nikkei'
+ ,'DJ_Indust'
+ ,'Goudprijs'
+ ,'Conjunct_Seizoenzuiver'
+ ,'Cons_vertrouw'
+ ,'Alg_consumptie_index_BE'
+ ,'Gem_rente_kasbon_1j')
+ ,1:51))
> y <- array(NA,dim=c(8,51),dimnames=list(c('BEL_20','Nikkei','DJ_Indust','Goudprijs','Conjunct_Seizoenzuiver','Cons_vertrouw','Alg_consumptie_index_BE','Gem_rente_kasbon_1j'),1:51))
> 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
BEL_20 Nikkei DJ_Indust Goudprijs Conjunct_Seizoenzuiver Cons_vertrouw
1 3484.74 13830.14 9349.44 7977 -5.6 6
2 3411.13 14153.22 9327.78 8241 -6.2 3
3 3288.18 15418.03 9753.63 8444 -7.1 2
4 3280.37 16666.97 10443.50 8490 -1.4 2
5 3173.95 16505.21 10853.87 8388 -0.1 2
6 3165.26 17135.96 10704.02 8099 -0.9 -8
7 3092.71 18033.25 11052.23 7984 0.0 0
8 3053.05 17671.00 10935.47 7786 0.1 -2
9 3181.96 17544.22 10714.03 8086 2.6 3
10 2999.93 17677.90 10394.48 9315 6.0 5
11 3249.57 18470.97 10817.90 9113 6.4 8
12 3210.52 18409.96 11251.20 9023 8.6 8
13 3030.29 18941.60 11281.26 9026 6.4 9
14 2803.47 19685.53 10539.68 9787 7.7 11
15 2767.63 19834.71 10483.39 9536 9.2 13
16 2882.60 19598.93 10947.43 9490 8.6 12
17 2863.36 17039.97 10580.27 9736 7.4 13
18 2897.06 16969.28 10582.92 9694 8.6 15
19 3012.61 16973.38 10654.41 9647 6.2 13
20 3142.95 16329.89 11014.51 9753 6.0 16
21 3032.93 16153.34 10967.87 10070 6.6 10
22 3045.78 15311.70 10433.56 10137 5.1 14
23 3110.52 14760.87 10665.78 9984 4.7 14
24 3013.24 14452.93 10666.71 9732 5.0 15
25 2987.10 13720.95 10682.74 9103 3.6 13
26 2995.55 13266.27 10777.22 9155 1.9 8
27 2833.18 12708.47 10052.60 9308 -0.1 7
28 2848.96 13411.84 10213.97 9394 -5.7 3
29 2794.83 13975.55 10546.82 9948 -5.6 3
30 2845.26 12974.89 10767.20 10177 -6.4 4
31 2915.02 12151.11 10444.50 10002 -7.7 4
32 2892.63 11576.21 10314.68 9728 -8.0 0
33 2604.42 9996.83 9042.56 10002 -11.9 -4
34 2641.65 10438.90 9220.75 10063 -15.4 -14
35 2659.81 10511.22 9721.84 10018 -15.5 -18
36 2638.53 10496.20 9978.53 9960 -13.4 -8
37 2720.25 10300.79 9923.81 10236 -10.9 -1
38 2745.88 9981.65 9892.56 10893 -10.8 1
39 2735.70 11448.79 10500.98 10756 -7.3 2
40 2811.70 11384.49 10179.35 10940 -6.5 0
41 2799.43 11717.46 10080.48 10997 -5.1 1
42 2555.28 10965.88 9492.44 10827 -5.3 0
43 2304.98 10352.27 8616.49 10166 -6.8 -1
44 2214.95 9751.20 8685.40 10186 -8.4 -3
45 2065.81 9354.01 8160.67 10457 -8.4 -3
46 1940.49 8792.50 8048.10 10368 -9.7 -3
47 2042.00 8721.14 8641.21 10244 -8.8 -4
48 1995.37 8692.94 8526.63 10511 -9.6 -8
49 1946.81 8570.73 8474.21 10812 -11.5 -9
50 1765.90 8538.47 7916.13 10738 -11.0 -13
51 1635.25 8169.75 7977.64 10171 -14.9 -18
Alg_consumptie_index_BE Gem_rente_kasbon_1j M1 M2 M3 M4 M5 M6 M7 M8 M9 M10
1 1.0 2.77 1 0 0 0 0 0 0 0 0 0
2 1.0 2.76 0 1 0 0 0 0 0 0 0 0
3 1.2 2.76 0 0 1 0 0 0 0 0 0 0
4 1.2 2.46 0 0 0 1 0 0 0 0 0 0
5 0.8 2.46 0 0 0 0 1 0 0 0 0 0
6 0.7 2.47 0 0 0 0 0 1 0 0 0 0
7 0.7 2.71 0 0 0 0 0 0 1 0 0 0
8 0.9 2.80 0 0 0 0 0 0 0 1 0 0
9 1.2 2.89 0 0 0 0 0 0 0 0 1 0
10 1.3 3.36 0 0 0 0 0 0 0 0 0 1
11 1.5 3.31 0 0 0 0 0 0 0 0 0 0
12 1.9 3.50 0 0 0 0 0 0 0 0 0 0
13 1.8 3.51 1 0 0 0 0 0 0 0 0 0
14 1.9 3.71 0 1 0 0 0 0 0 0 0 0
15 2.2 3.71 0 0 1 0 0 0 0 0 0 0
16 2.1 3.71 0 0 0 1 0 0 0 0 0 0
17 2.2 4.21 0 0 0 0 1 0 0 0 0 0
18 2.7 4.21 0 0 0 0 0 1 0 0 0 0
19 2.8 4.21 0 0 0 0 0 0 1 0 0 0
20 2.9 4.50 0 0 0 0 0 0 0 1 0 0
21 3.4 4.51 0 0 0 0 0 0 0 0 1 0
22 3.0 4.51 0 0 0 0 0 0 0 0 0 1
23 3.1 4.51 0 0 0 0 0 0 0 0 0 0
24 2.5 4.32 0 0 0 0 0 0 0 0 0 0
25 2.2 4.02 1 0 0 0 0 0 0 0 0 0
26 2.3 4.02 0 1 0 0 0 0 0 0 0 0
27 2.1 3.85 0 0 1 0 0 0 0 0 0 0
28 2.8 3.84 0 0 0 1 0 0 0 0 0 0
29 3.1 4.02 0 0 0 0 1 0 0 0 0 0
30 2.9 3.82 0 0 0 0 0 1 0 0 0 0
31 2.6 3.75 0 0 0 0 0 0 1 0 0 0
32 2.7 3.74 0 0 0 0 0 0 0 1 0 0
33 2.3 3.14 0 0 0 0 0 0 0 0 1 0
34 2.3 2.91 0 0 0 0 0 0 0 0 0 1
35 2.1 2.84 0 0 0 0 0 0 0 0 0 0
36 2.2 2.85 0 0 0 0 0 0 0 0 0 0
37 2.9 2.85 1 0 0 0 0 0 0 0 0 0
38 2.6 3.08 0 1 0 0 0 0 0 0 0 0
39 2.7 3.30 0 0 1 0 0 0 0 0 0 0
40 1.8 3.29 0 0 0 1 0 0 0 0 0 0
41 1.3 3.26 0 0 0 0 1 0 0 0 0 0
42 0.9 3.26 0 0 0 0 0 1 0 0 0 0
43 1.3 3.11 0 0 0 0 0 0 1 0 0 0
44 1.3 2.84 0 0 0 0 0 0 0 1 0 0
45 1.3 2.71 0 0 0 0 0 0 0 0 1 0
46 1.3 2.69 0 0 0 0 0 0 0 0 0 1
47 1.1 2.65 0 0 0 0 0 0 0 0 0 0
48 1.4 2.57 0 0 0 0 0 0 0 0 0 0
49 1.2 2.32 1 0 0 0 0 0 0 0 0 0
50 1.7 2.12 0 1 0 0 0 0 0 0 0 0
51 1.8 2.05 0 0 1 0 0 0 0 0 0 0
M11
1 0
2 0
3 0
4 0
5 0
6 0
7 0
8 0
9 0
10 0
11 1
12 0
13 0
14 0
15 0
16 0
17 0
18 0
19 0
20 0
21 0
22 0
23 1
24 0
25 0
26 0
27 0
28 0
29 0
30 0
31 0
32 0
33 0
34 0
35 1
36 0
37 0
38 0
39 0
40 0
41 0
42 0
43 0
44 0
45 0
46 0
47 1
48 0
49 0
50 0
51 0
> k <- length(x[1,])
> df <- as.data.frame(x)
> (mylm <- lm(df))
Call:
lm(formula = df)
Coefficients:
(Intercept) Nikkei DJ_Indust
688.09880 0.04622 0.27803
Goudprijs Conjunct_Seizoenzuiver Cons_vertrouw
-0.17878 -57.55484 33.67602
Alg_consumptie_index_BE Gem_rente_kasbon_1j M1
-131.10901 131.83865 -34.42476
M2 M3 M4
56.98807 -72.54450 -24.12459
M5 M6 M7
-68.21470 -14.56166 -131.28579
M8 M9 M10
-116.52423 65.12671 96.69511
M11
87.69527
> (mysum <- summary(mylm))
Call:
lm(formula = df)
Residuals:
Min 1Q Median 3Q Max
-284.92 -100.64 16.01 93.64 260.10
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 688.09880 899.14335 0.765 0.449711
Nikkei 0.04622 0.02709 1.706 0.097682 .
DJ_Indust 0.27803 0.06205 4.481 8.92e-05 ***
Goudprijs -0.17878 0.05149 -3.472 0.001502 **
Conjunct_Seizoenzuiver -57.55484 12.58878 -4.572 6.86e-05 ***
Cons_vertrouw 33.67602 8.68637 3.877 0.000494 ***
Alg_consumptie_index_BE -131.10901 75.27088 -1.742 0.091145 .
Gem_rente_kasbon_1j 131.83865 130.18537 1.013 0.318800
M1 -34.42476 119.64691 -0.288 0.775417
M2 56.98807 118.58659 0.481 0.634098
M3 -72.54450 120.30205 -0.603 0.550746
M4 -24.12459 126.55526 -0.191 0.850024
M5 -68.21470 127.07337 -0.537 0.595109
M6 -14.56166 124.99330 -0.116 0.907985
M7 -131.28579 126.34021 -1.039 0.306526
M8 -116.52423 124.03408 -0.939 0.354536
M9 65.12671 124.00618 0.525 0.603072
M10 96.69511 129.14292 0.749 0.459479
M11 87.69527 122.46031 0.716 0.479116
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 169.6 on 32 degrees of freedom
Multiple R-squared: 0.901, Adjusted R-squared: 0.8453
F-statistic: 16.18 on 18 and 32 DF, p-value: 2.856e-11
> 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.1886796 0.3773591 0.8113204
[2,] 0.3476281 0.6952562 0.6523719
[3,] 0.4532082 0.9064164 0.5467918
[4,] 0.3185107 0.6370215 0.6814893
[5,] 0.4299946 0.8599891 0.5700054
[6,] 0.2827890 0.5655780 0.7172110
[7,] 0.3742126 0.7484253 0.6257874
[8,] 0.2669282 0.5338565 0.7330718
> postscript(file="/var/www/html/rcomp/tmp/1ya0q1291646817.ps",horizontal=F,onefile=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/2rjhb1291646817.ps",horizontal=F,onefile=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/3rjhb1291646817.ps",horizontal=F,onefile=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/4rjhb1291646817.ps",horizontal=F,onefile=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/51azw1291646817.ps",horizontal=F,onefile=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 = 51
Frequency = 1
1 2 3 4 5 6
260.0991192 201.1772153 75.2939088 145.3736179 -19.4330309 155.3550827
7 8 9 10 11 12
-208.5642687 -161.7146748 -90.4215002 77.8453720 100.8019380 169.7212990
13 14 15 16 17 18
-183.2024997 -199.3741593 -83.4863745 -157.2478019 -23.6169530 18.7199655
19 20 21 22 23 24
164.8569812 91.3465842 178.3005002 85.5336030 82.9034299 -27.7851333
25 26 27 28 29 30
-115.5836268 -110.8561319 25.6633519 -163.4841097 -171.7198324 -228.6005751
31 32 33 34 35 36
-50.5339739 57.8609182 0.5687336 112.8031474 101.2127492 -117.5170894
37 38 39 40 41 42
72.1468615 16.0073553 25.7721435 175.3582937 214.7698163 54.5255270
43 44 45 46 47 48
94.2412614 12.5071724 -88.4477336 -276.1821224 -284.9181171 -24.4190763
49 50 51
-33.4598542 93.0457207 -43.2430297
> postscript(file="/var/www/html/rcomp/tmp/61azw1291646817.ps",horizontal=F,onefile=F,pagecentre=F,paper="special",width=8.3333333333333,height=5.5555555555556)
> dum <- cbind(lag(myerror,k=1),myerror)
> dum
Time Series:
Start = 0
End = 51
Frequency = 1
lag(myerror, k = 1) myerror
0 260.0991192 NA
1 201.1772153 260.0991192
2 75.2939088 201.1772153
3 145.3736179 75.2939088
4 -19.4330309 145.3736179
5 155.3550827 -19.4330309
6 -208.5642687 155.3550827
7 -161.7146748 -208.5642687
8 -90.4215002 -161.7146748
9 77.8453720 -90.4215002
10 100.8019380 77.8453720
11 169.7212990 100.8019380
12 -183.2024997 169.7212990
13 -199.3741593 -183.2024997
14 -83.4863745 -199.3741593
15 -157.2478019 -83.4863745
16 -23.6169530 -157.2478019
17 18.7199655 -23.6169530
18 164.8569812 18.7199655
19 91.3465842 164.8569812
20 178.3005002 91.3465842
21 85.5336030 178.3005002
22 82.9034299 85.5336030
23 -27.7851333 82.9034299
24 -115.5836268 -27.7851333
25 -110.8561319 -115.5836268
26 25.6633519 -110.8561319
27 -163.4841097 25.6633519
28 -171.7198324 -163.4841097
29 -228.6005751 -171.7198324
30 -50.5339739 -228.6005751
31 57.8609182 -50.5339739
32 0.5687336 57.8609182
33 112.8031474 0.5687336
34 101.2127492 112.8031474
35 -117.5170894 101.2127492
36 72.1468615 -117.5170894
37 16.0073553 72.1468615
38 25.7721435 16.0073553
39 175.3582937 25.7721435
40 214.7698163 175.3582937
41 54.5255270 214.7698163
42 94.2412614 54.5255270
43 12.5071724 94.2412614
44 -88.4477336 12.5071724
45 -276.1821224 -88.4477336
46 -284.9181171 -276.1821224
47 -24.4190763 -284.9181171
48 -33.4598542 -24.4190763
49 93.0457207 -33.4598542
50 -43.2430297 93.0457207
51 NA -43.2430297
> dum1 <- dum[2:length(myerror),]
> dum1
lag(myerror, k = 1) myerror
[1,] 201.1772153 260.0991192
[2,] 75.2939088 201.1772153
[3,] 145.3736179 75.2939088
[4,] -19.4330309 145.3736179
[5,] 155.3550827 -19.4330309
[6,] -208.5642687 155.3550827
[7,] -161.7146748 -208.5642687
[8,] -90.4215002 -161.7146748
[9,] 77.8453720 -90.4215002
[10,] 100.8019380 77.8453720
[11,] 169.7212990 100.8019380
[12,] -183.2024997 169.7212990
[13,] -199.3741593 -183.2024997
[14,] -83.4863745 -199.3741593
[15,] -157.2478019 -83.4863745
[16,] -23.6169530 -157.2478019
[17,] 18.7199655 -23.6169530
[18,] 164.8569812 18.7199655
[19,] 91.3465842 164.8569812
[20,] 178.3005002 91.3465842
[21,] 85.5336030 178.3005002
[22,] 82.9034299 85.5336030
[23,] -27.7851333 82.9034299
[24,] -115.5836268 -27.7851333
[25,] -110.8561319 -115.5836268
[26,] 25.6633519 -110.8561319
[27,] -163.4841097 25.6633519
[28,] -171.7198324 -163.4841097
[29,] -228.6005751 -171.7198324
[30,] -50.5339739 -228.6005751
[31,] 57.8609182 -50.5339739
[32,] 0.5687336 57.8609182
[33,] 112.8031474 0.5687336
[34,] 101.2127492 112.8031474
[35,] -117.5170894 101.2127492
[36,] 72.1468615 -117.5170894
[37,] 16.0073553 72.1468615
[38,] 25.7721435 16.0073553
[39,] 175.3582937 25.7721435
[40,] 214.7698163 175.3582937
[41,] 54.5255270 214.7698163
[42,] 94.2412614 54.5255270
[43,] 12.5071724 94.2412614
[44,] -88.4477336 12.5071724
[45,] -276.1821224 -88.4477336
[46,] -284.9181171 -276.1821224
[47,] -24.4190763 -284.9181171
[48,] -33.4598542 -24.4190763
[49,] 93.0457207 -33.4598542
[50,] -43.2430297 93.0457207
> z <- as.data.frame(dum1)
> z
lag(myerror, k = 1) myerror
1 201.1772153 260.0991192
2 75.2939088 201.1772153
3 145.3736179 75.2939088
4 -19.4330309 145.3736179
5 155.3550827 -19.4330309
6 -208.5642687 155.3550827
7 -161.7146748 -208.5642687
8 -90.4215002 -161.7146748
9 77.8453720 -90.4215002
10 100.8019380 77.8453720
11 169.7212990 100.8019380
12 -183.2024997 169.7212990
13 -199.3741593 -183.2024997
14 -83.4863745 -199.3741593
15 -157.2478019 -83.4863745
16 -23.6169530 -157.2478019
17 18.7199655 -23.6169530
18 164.8569812 18.7199655
19 91.3465842 164.8569812
20 178.3005002 91.3465842
21 85.5336030 178.3005002
22 82.9034299 85.5336030
23 -27.7851333 82.9034299
24 -115.5836268 -27.7851333
25 -110.8561319 -115.5836268
26 25.6633519 -110.8561319
27 -163.4841097 25.6633519
28 -171.7198324 -163.4841097
29 -228.6005751 -171.7198324
30 -50.5339739 -228.6005751
31 57.8609182 -50.5339739
32 0.5687336 57.8609182
33 112.8031474 0.5687336
34 101.2127492 112.8031474
35 -117.5170894 101.2127492
36 72.1468615 -117.5170894
37 16.0073553 72.1468615
38 25.7721435 16.0073553
39 175.3582937 25.7721435
40 214.7698163 175.3582937
41 54.5255270 214.7698163
42 94.2412614 54.5255270
43 12.5071724 94.2412614
44 -88.4477336 12.5071724
45 -276.1821224 -88.4477336
46 -284.9181171 -276.1821224
47 -24.4190763 -284.9181171
48 -33.4598542 -24.4190763
49 93.0457207 -33.4598542
50 -43.2430297 93.0457207
> 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/7c1yz1291646817.ps",horizontal=F,onefile=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/8c1yz1291646817.ps",horizontal=F,onefile=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/9ntx21291646817.ps",horizontal=F,onefile=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/10ntx21291646817.ps",horizontal=F,onefile=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/11qteq1291646817.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/12cucw1291646817.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/130v971291646817.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/14t48s1291646817.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/15e46y1291646817.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/16awm71291646817.tab")
+ }
>
> try(system("convert tmp/1ya0q1291646817.ps tmp/1ya0q1291646817.png",intern=TRUE))
character(0)
> try(system("convert tmp/2rjhb1291646817.ps tmp/2rjhb1291646817.png",intern=TRUE))
character(0)
> try(system("convert tmp/3rjhb1291646817.ps tmp/3rjhb1291646817.png",intern=TRUE))
character(0)
> try(system("convert tmp/4rjhb1291646817.ps tmp/4rjhb1291646817.png",intern=TRUE))
character(0)
> try(system("convert tmp/51azw1291646817.ps tmp/51azw1291646817.png",intern=TRUE))
character(0)
> try(system("convert tmp/61azw1291646817.ps tmp/61azw1291646817.png",intern=TRUE))
character(0)
> try(system("convert tmp/7c1yz1291646817.ps tmp/7c1yz1291646817.png",intern=TRUE))
character(0)
> try(system("convert tmp/8c1yz1291646817.ps tmp/8c1yz1291646817.png",intern=TRUE))
character(0)
> try(system("convert tmp/9ntx21291646817.ps tmp/9ntx21291646817.png",intern=TRUE))
character(0)
> try(system("convert tmp/10ntx21291646817.ps tmp/10ntx21291646817.png",intern=TRUE))
character(0)
>
>
> proc.time()
user system elapsed
2.375 1.685 6.192