R version 2.12.0 (2010-10-15)
Copyright (C) 2010 The R Foundation for Statistical Computing
ISBN 3-900051-07-0
Platform: i486-pc-linux-gnu (32-bit)
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 = 'Do not include Seasonal 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
> 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
1 1.0 2.77
2 1.0 2.76
3 1.2 2.76
4 1.2 2.46
5 0.8 2.46
6 0.7 2.47
7 0.7 2.71
8 0.9 2.80
9 1.2 2.89
10 1.3 3.36
11 1.5 3.31
12 1.9 3.50
13 1.8 3.51
14 1.9 3.71
15 2.2 3.71
16 2.1 3.71
17 2.2 4.21
18 2.7 4.21
19 2.8 4.21
20 2.9 4.50
21 3.4 4.51
22 3.0 4.51
23 3.1 4.51
24 2.5 4.32
25 2.2 4.02
26 2.3 4.02
27 2.1 3.85
28 2.8 3.84
29 3.1 4.02
30 2.9 3.82
31 2.6 3.75
32 2.7 3.74
33 2.3 3.14
34 2.3 2.91
35 2.1 2.84
36 2.2 2.85
37 2.9 2.85
38 2.6 3.08
39 2.7 3.30
40 1.8 3.29
41 1.3 3.26
42 0.9 3.26
43 1.3 3.11
44 1.3 2.84
45 1.3 2.71
46 1.3 2.69
47 1.1 2.65
48 1.4 2.57
49 1.2 2.32
50 1.7 2.12
51 1.8 2.05
> k <- length(x[1,])
> df <- as.data.frame(x)
> (mylm <- lm(df))
Call:
lm(formula = df)
Coefficients:
(Intercept) Nikkei DJ_Indust
961.21319 0.04524 0.24917
Goudprijs Conjunct_Seizoenzuiver Cons_vertrouw
-0.17784 -51.37497 29.90473
Alg_consumptie_index_BE Gem_rente_kasbon_1j
-105.78782 127.07989
> (mysum <- summary(mylm))
Call:
lm(formula = df)
Residuals:
Min 1Q Median 3Q Max
-288.57 -127.94 -15.63 155.39 297.90
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 961.21319 756.28240 1.271 0.210572
Nikkei 0.04524 0.02374 1.905 0.063444 .
DJ_Indust 0.24917 0.05338 4.667 2.98e-05 ***
Goudprijs -0.17784 0.04598 -3.868 0.000367 ***
Conjunct_Seizoenzuiver -51.37497 11.30025 -4.546 4.40e-05 ***
Cons_vertrouw 29.90473 7.37313 4.056 0.000206 ***
Alg_consumptie_index_BE -105.78782 63.50633 -1.666 0.103026
Gem_rente_kasbon_1j 127.07989 106.53369 1.193 0.239465
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 164.3 on 43 degrees of freedom
Multiple R-squared: 0.8752, Adjusted R-squared: 0.8548
F-statistic: 43.07 on 7 and 43 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.6047610 0.790477950 0.395238975
[2,] 0.5447600 0.910480044 0.455240022
[3,] 0.4655854 0.931170874 0.534414563
[4,] 0.5255223 0.948955369 0.474477684
[5,] 0.5480252 0.903949671 0.451974836
[6,] 0.4483780 0.896756080 0.551621960
[7,] 0.6616614 0.676677215 0.338338607
[8,] 0.6492819 0.701436137 0.350718069
[9,] 0.5780208 0.843958400 0.421979200
[10,] 0.5375925 0.924815082 0.462407541
[11,] 0.4380995 0.876199023 0.561900488
[12,] 0.3626250 0.725250057 0.637374971
[13,] 0.3154861 0.630972229 0.684513886
[14,] 0.3402580 0.680516076 0.659741962
[15,] 0.4415999 0.883199769 0.558400116
[16,] 0.4443424 0.888684790 0.555657605
[17,] 0.7509174 0.498165193 0.249082596
[18,] 0.8374760 0.325048083 0.162524041
[19,] 0.8475973 0.304805499 0.152402750
[20,] 0.9505617 0.098876578 0.049438289
[21,] 0.9762218 0.047556332 0.023778166
[22,] 0.9624144 0.075171246 0.037585623
[23,] 0.9597755 0.080448990 0.040224495
[24,] 0.9342709 0.131458188 0.065729094
[25,] 0.9473586 0.105282733 0.052641367
[26,] 0.9108816 0.178236837 0.089118419
[27,] 0.9445111 0.110977737 0.055488869
[28,] 0.9898297 0.020340566 0.010170283
[29,] 0.9986712 0.002657589 0.001328794
[30,] 0.9930785 0.013843078 0.006921539
> postscript(file="/var/www/rcomp/tmp/1lnvs1291646534.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/rcomp/tmp/2lnvs1291646534.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/rcomp/tmp/3lnvs1291646534.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/rcomp/tmp/4eevd1291646534.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/rcomp/tmp/5eevd1291646534.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
273.6144200 297.8969951 52.5498455 155.4901716 -39.5314124 155.2846108
7 8 9 10 11 12
-288.5702192 -243.2959662 -0.9090500 174.9215593 205.6058291 176.5398941
13 14 15 16 17 18
-189.4753873 -137.6901760 -161.9021486 -171.5703616 -84.3335014 -0.8310211
19 20 21 22 23 24
35.4501136 -2.2384724 225.6020204 182.5772491 177.1917143 -5.0261063
25 26 27 28 29 30
-119.6392880 -32.1497841 -33.9233796 -167.6339016 -217.6754152 -202.9110237
31 32 33 34 35 36
-136.2300892 -32.9395507 69.1870366 201.3316302 185.5806591 -91.1457705
37 38 39 40 41 42
55.2889448 104.3485543 -15.6346685 183.1085205 183.4872020 66.9385276
43 44 45 46 47 48
-40.6794890 -105.2114422 -40.9217383 -192.8659518 -197.8970649 -46.7958368
49 50 51
-80.3305613 89.7307472 -171.7674673
> postscript(file="/var/www/rcomp/tmp/6eevd1291646534.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 273.6144200 NA
1 297.8969951 273.6144200
2 52.5498455 297.8969951
3 155.4901716 52.5498455
4 -39.5314124 155.4901716
5 155.2846108 -39.5314124
6 -288.5702192 155.2846108
7 -243.2959662 -288.5702192
8 -0.9090500 -243.2959662
9 174.9215593 -0.9090500
10 205.6058291 174.9215593
11 176.5398941 205.6058291
12 -189.4753873 176.5398941
13 -137.6901760 -189.4753873
14 -161.9021486 -137.6901760
15 -171.5703616 -161.9021486
16 -84.3335014 -171.5703616
17 -0.8310211 -84.3335014
18 35.4501136 -0.8310211
19 -2.2384724 35.4501136
20 225.6020204 -2.2384724
21 182.5772491 225.6020204
22 177.1917143 182.5772491
23 -5.0261063 177.1917143
24 -119.6392880 -5.0261063
25 -32.1497841 -119.6392880
26 -33.9233796 -32.1497841
27 -167.6339016 -33.9233796
28 -217.6754152 -167.6339016
29 -202.9110237 -217.6754152
30 -136.2300892 -202.9110237
31 -32.9395507 -136.2300892
32 69.1870366 -32.9395507
33 201.3316302 69.1870366
34 185.5806591 201.3316302
35 -91.1457705 185.5806591
36 55.2889448 -91.1457705
37 104.3485543 55.2889448
38 -15.6346685 104.3485543
39 183.1085205 -15.6346685
40 183.4872020 183.1085205
41 66.9385276 183.4872020
42 -40.6794890 66.9385276
43 -105.2114422 -40.6794890
44 -40.9217383 -105.2114422
45 -192.8659518 -40.9217383
46 -197.8970649 -192.8659518
47 -46.7958368 -197.8970649
48 -80.3305613 -46.7958368
49 89.7307472 -80.3305613
50 -171.7674673 89.7307472
51 NA -171.7674673
> dum1 <- dum[2:length(myerror),]
> dum1
lag(myerror, k = 1) myerror
[1,] 297.8969951 273.6144200
[2,] 52.5498455 297.8969951
[3,] 155.4901716 52.5498455
[4,] -39.5314124 155.4901716
[5,] 155.2846108 -39.5314124
[6,] -288.5702192 155.2846108
[7,] -243.2959662 -288.5702192
[8,] -0.9090500 -243.2959662
[9,] 174.9215593 -0.9090500
[10,] 205.6058291 174.9215593
[11,] 176.5398941 205.6058291
[12,] -189.4753873 176.5398941
[13,] -137.6901760 -189.4753873
[14,] -161.9021486 -137.6901760
[15,] -171.5703616 -161.9021486
[16,] -84.3335014 -171.5703616
[17,] -0.8310211 -84.3335014
[18,] 35.4501136 -0.8310211
[19,] -2.2384724 35.4501136
[20,] 225.6020204 -2.2384724
[21,] 182.5772491 225.6020204
[22,] 177.1917143 182.5772491
[23,] -5.0261063 177.1917143
[24,] -119.6392880 -5.0261063
[25,] -32.1497841 -119.6392880
[26,] -33.9233796 -32.1497841
[27,] -167.6339016 -33.9233796
[28,] -217.6754152 -167.6339016
[29,] -202.9110237 -217.6754152
[30,] -136.2300892 -202.9110237
[31,] -32.9395507 -136.2300892
[32,] 69.1870366 -32.9395507
[33,] 201.3316302 69.1870366
[34,] 185.5806591 201.3316302
[35,] -91.1457705 185.5806591
[36,] 55.2889448 -91.1457705
[37,] 104.3485543 55.2889448
[38,] -15.6346685 104.3485543
[39,] 183.1085205 -15.6346685
[40,] 183.4872020 183.1085205
[41,] 66.9385276 183.4872020
[42,] -40.6794890 66.9385276
[43,] -105.2114422 -40.6794890
[44,] -40.9217383 -105.2114422
[45,] -192.8659518 -40.9217383
[46,] -197.8970649 -192.8659518
[47,] -46.7958368 -197.8970649
[48,] -80.3305613 -46.7958368
[49,] 89.7307472 -80.3305613
[50,] -171.7674673 89.7307472
> z <- as.data.frame(dum1)
> z
lag(myerror, k = 1) myerror
1 297.8969951 273.6144200
2 52.5498455 297.8969951
3 155.4901716 52.5498455
4 -39.5314124 155.4901716
5 155.2846108 -39.5314124
6 -288.5702192 155.2846108
7 -243.2959662 -288.5702192
8 -0.9090500 -243.2959662
9 174.9215593 -0.9090500
10 205.6058291 174.9215593
11 176.5398941 205.6058291
12 -189.4753873 176.5398941
13 -137.6901760 -189.4753873
14 -161.9021486 -137.6901760
15 -171.5703616 -161.9021486
16 -84.3335014 -171.5703616
17 -0.8310211 -84.3335014
18 35.4501136 -0.8310211
19 -2.2384724 35.4501136
20 225.6020204 -2.2384724
21 182.5772491 225.6020204
22 177.1917143 182.5772491
23 -5.0261063 177.1917143
24 -119.6392880 -5.0261063
25 -32.1497841 -119.6392880
26 -33.9233796 -32.1497841
27 -167.6339016 -33.9233796
28 -217.6754152 -167.6339016
29 -202.9110237 -217.6754152
30 -136.2300892 -202.9110237
31 -32.9395507 -136.2300892
32 69.1870366 -32.9395507
33 201.3316302 69.1870366
34 185.5806591 201.3316302
35 -91.1457705 185.5806591
36 55.2889448 -91.1457705
37 104.3485543 55.2889448
38 -15.6346685 104.3485543
39 183.1085205 -15.6346685
40 183.4872020 183.1085205
41 66.9385276 183.4872020
42 -40.6794890 66.9385276
43 -105.2114422 -40.6794890
44 -40.9217383 -105.2114422
45 -192.8659518 -40.9217383
46 -197.8970649 -192.8659518
47 -46.7958368 -197.8970649
48 -80.3305613 -46.7958368
49 89.7307472 -80.3305613
50 -171.7674673 89.7307472
> 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/rcomp/tmp/7p5cf1291646534.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/rcomp/tmp/8zfbi1291646534.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/rcomp/tmp/9zfbi1291646534.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/rcomp/tmp/10zfbi1291646534.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/rcomp/createtable file can be downloaded at http://www.wessa.net/cretab
> load(file="/var/www/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/rcomp/tmp/11lxs61291646534.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/rcomp/tmp/126y8c1291646534.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/rcomp/tmp/13vzn61291646534.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/rcomp/tmp/146qm91291646534.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/rcomp/tmp/15rq3x1291646534.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/rcomp/tmp/1650j51291646534.tab")
+ }
>
> try(system("convert tmp/1lnvs1291646534.ps tmp/1lnvs1291646534.png",intern=TRUE))
character(0)
> try(system("convert tmp/2lnvs1291646534.ps tmp/2lnvs1291646534.png",intern=TRUE))
character(0)
> try(system("convert tmp/3lnvs1291646534.ps tmp/3lnvs1291646534.png",intern=TRUE))
character(0)
> try(system("convert tmp/4eevd1291646534.ps tmp/4eevd1291646534.png",intern=TRUE))
character(0)
> try(system("convert tmp/5eevd1291646534.ps tmp/5eevd1291646534.png",intern=TRUE))
character(0)
> try(system("convert tmp/6eevd1291646534.ps tmp/6eevd1291646534.png",intern=TRUE))
character(0)
> try(system("convert tmp/7p5cf1291646534.ps tmp/7p5cf1291646534.png",intern=TRUE))
character(0)
> try(system("convert tmp/8zfbi1291646534.ps tmp/8zfbi1291646534.png",intern=TRUE))
character(0)
> try(system("convert tmp/9zfbi1291646534.ps tmp/9zfbi1291646534.png",intern=TRUE))
character(0)
> try(system("convert tmp/10zfbi1291646534.ps tmp/10zfbi1291646534.png",intern=TRUE))
character(0)
>
>
> proc.time()
user system elapsed
3.000 1.760 4.733