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 = '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
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 t
1 1.0 2.77 1
2 1.0 2.76 2
3 1.2 2.76 3
4 1.2 2.46 4
5 0.8 2.46 5
6 0.7 2.47 6
7 0.7 2.71 7
8 0.9 2.80 8
9 1.2 2.89 9
10 1.3 3.36 10
11 1.5 3.31 11
12 1.9 3.50 12
13 1.8 3.51 13
14 1.9 3.71 14
15 2.2 3.71 15
16 2.1 3.71 16
17 2.2 4.21 17
18 2.7 4.21 18
19 2.8 4.21 19
20 2.9 4.50 20
21 3.4 4.51 21
22 3.0 4.51 22
23 3.1 4.51 23
24 2.5 4.32 24
25 2.2 4.02 25
26 2.3 4.02 26
27 2.1 3.85 27
28 2.8 3.84 28
29 3.1 4.02 29
30 2.9 3.82 30
31 2.6 3.75 31
32 2.7 3.74 32
33 2.3 3.14 33
34 2.3 2.91 34
35 2.1 2.84 35
36 2.2 2.85 36
37 2.9 2.85 37
38 2.6 3.08 38
39 2.7 3.30 39
40 1.8 3.29 40
41 1.3 3.26 41
42 0.9 3.26 42
43 1.3 3.11 43
44 1.3 2.84 44
45 1.3 2.71 45
46 1.3 2.69 46
47 1.1 2.65 47
48 1.4 2.57 48
49 1.2 2.32 49
50 1.7 2.12 50
51 1.8 2.05 51
> k <- length(x[1,])
> df <- as.data.frame(x)
> (mylm <- lm(df))
Call:
lm(formula = df)
Coefficients:
(Intercept) Nikkei DJ_Indust
623.18928 -0.08298 0.26584
Goudprijs Conjunct_Seizoenzuiver Cons_vertrouw
0.13217 -13.12259 6.30505
Alg_consumptie_index_BE Gem_rente_kasbon_1j t
-96.56682 146.91292 -38.08906
> (mysum <- summary(mylm))
Call:
lm(formula = df)
Residuals:
Min 1Q Median 3Q Max
-188.639 -45.650 2.374 41.090 189.984
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 623.18928 412.10320 1.512 0.13797
Nikkei -0.08298 0.01802 -4.605 3.80e-05 ***
DJ_Indust 0.26584 0.02904 9.154 1.47e-11 ***
Goudprijs 0.13217 0.03937 3.357 0.00168 **
Conjunct_Seizoenzuiver -13.12259 7.19525 -1.824 0.07531 .
Cons_vertrouw 6.30505 4.62649 1.363 0.18020
Alg_consumptie_index_BE -96.56682 34.50452 -2.799 0.00772 **
Gem_rente_kasbon_1j 146.91292 57.89514 2.538 0.01496 *
t -38.08906 3.73919 -10.186 6.44e-13 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 89.24 on 42 degrees of freedom
Multiple R-squared: 0.964, Adjusted R-squared: 0.9572
F-statistic: 140.7 on 8 and 42 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.8718270 0.25634592 0.12817296
[2,] 0.7953575 0.40928491 0.20464246
[3,] 0.8422984 0.31540316 0.15770158
[4,] 0.8496177 0.30076452 0.15038226
[5,] 0.7856298 0.42874046 0.21437023
[6,] 0.8443684 0.31126321 0.15563161
[7,] 0.8213821 0.35723584 0.17861792
[8,] 0.8581058 0.28378845 0.14189423
[9,] 0.8561976 0.28760477 0.14380239
[10,] 0.8091777 0.38164452 0.19082226
[11,] 0.8136478 0.37270441 0.18635220
[12,] 0.7701667 0.45966651 0.22983325
[13,] 0.7051387 0.58972256 0.29486128
[14,] 0.6172830 0.76543390 0.38271695
[15,] 0.5186627 0.96267468 0.48133734
[16,] 0.4615735 0.92314690 0.53842655
[17,] 0.3782689 0.75653773 0.62173114
[18,] 0.3371191 0.67423815 0.66288093
[19,] 0.8606062 0.27878759 0.13939379
[20,] 0.9576100 0.08477994 0.04238997
[21,] 0.9354778 0.12904436 0.06452218
[22,] 0.9246705 0.15065910 0.07532955
[23,] 0.8860690 0.22786200 0.11393100
[24,] 0.8183839 0.36323227 0.18161613
[25,] 0.9386060 0.12278804 0.06139402
[26,] 0.8838605 0.23227904 0.11613952
[27,] 0.8462975 0.30740504 0.15370252
[28,] 0.9589741 0.08205173 0.04102587
> postscript(file="/var/www/html/rcomp/tmp/1wo091291646919.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/2wo091291646919.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/3wo091291646919.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/46xit1291646919.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/56xit1291646919.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
85.856374 60.521959 -45.609982 17.709353 -181.223659 19.977176
7 8 9 10 11 12
-91.280613 -45.689998 147.034075 -90.686608 189.984027 120.250250
13 14 15 16 17 18
-32.461252 -78.174363 20.636723 25.624322 -188.639454 -66.449475
19 20 21 22 23 24
65.509498 16.312724 -7.249921 23.500947 43.516282 -40.564182
25 26 27 28 29 30
-1.142250 -5.450163 -17.887298 60.884764 -66.245814 -156.351938
31 32 33 34 35 36
-43.692932 27.433545 -28.222446 79.260493 29.127606 -43.187717
37 38 39 40 41 42
94.741286 -20.608911 2.374403 109.977472 150.368148 25.784874
43 44 45 46 47 48
130.164972 38.662888 17.426155 -88.833624 -91.760744 -52.026022
49 50 51
-99.699162 22.390478 -41.892266
> postscript(file="/var/www/html/rcomp/tmp/66xit1291646919.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 85.856374 NA
1 60.521959 85.856374
2 -45.609982 60.521959
3 17.709353 -45.609982
4 -181.223659 17.709353
5 19.977176 -181.223659
6 -91.280613 19.977176
7 -45.689998 -91.280613
8 147.034075 -45.689998
9 -90.686608 147.034075
10 189.984027 -90.686608
11 120.250250 189.984027
12 -32.461252 120.250250
13 -78.174363 -32.461252
14 20.636723 -78.174363
15 25.624322 20.636723
16 -188.639454 25.624322
17 -66.449475 -188.639454
18 65.509498 -66.449475
19 16.312724 65.509498
20 -7.249921 16.312724
21 23.500947 -7.249921
22 43.516282 23.500947
23 -40.564182 43.516282
24 -1.142250 -40.564182
25 -5.450163 -1.142250
26 -17.887298 -5.450163
27 60.884764 -17.887298
28 -66.245814 60.884764
29 -156.351938 -66.245814
30 -43.692932 -156.351938
31 27.433545 -43.692932
32 -28.222446 27.433545
33 79.260493 -28.222446
34 29.127606 79.260493
35 -43.187717 29.127606
36 94.741286 -43.187717
37 -20.608911 94.741286
38 2.374403 -20.608911
39 109.977472 2.374403
40 150.368148 109.977472
41 25.784874 150.368148
42 130.164972 25.784874
43 38.662888 130.164972
44 17.426155 38.662888
45 -88.833624 17.426155
46 -91.760744 -88.833624
47 -52.026022 -91.760744
48 -99.699162 -52.026022
49 22.390478 -99.699162
50 -41.892266 22.390478
51 NA -41.892266
> dum1 <- dum[2:length(myerror),]
> dum1
lag(myerror, k = 1) myerror
[1,] 60.521959 85.856374
[2,] -45.609982 60.521959
[3,] 17.709353 -45.609982
[4,] -181.223659 17.709353
[5,] 19.977176 -181.223659
[6,] -91.280613 19.977176
[7,] -45.689998 -91.280613
[8,] 147.034075 -45.689998
[9,] -90.686608 147.034075
[10,] 189.984027 -90.686608
[11,] 120.250250 189.984027
[12,] -32.461252 120.250250
[13,] -78.174363 -32.461252
[14,] 20.636723 -78.174363
[15,] 25.624322 20.636723
[16,] -188.639454 25.624322
[17,] -66.449475 -188.639454
[18,] 65.509498 -66.449475
[19,] 16.312724 65.509498
[20,] -7.249921 16.312724
[21,] 23.500947 -7.249921
[22,] 43.516282 23.500947
[23,] -40.564182 43.516282
[24,] -1.142250 -40.564182
[25,] -5.450163 -1.142250
[26,] -17.887298 -5.450163
[27,] 60.884764 -17.887298
[28,] -66.245814 60.884764
[29,] -156.351938 -66.245814
[30,] -43.692932 -156.351938
[31,] 27.433545 -43.692932
[32,] -28.222446 27.433545
[33,] 79.260493 -28.222446
[34,] 29.127606 79.260493
[35,] -43.187717 29.127606
[36,] 94.741286 -43.187717
[37,] -20.608911 94.741286
[38,] 2.374403 -20.608911
[39,] 109.977472 2.374403
[40,] 150.368148 109.977472
[41,] 25.784874 150.368148
[42,] 130.164972 25.784874
[43,] 38.662888 130.164972
[44,] 17.426155 38.662888
[45,] -88.833624 17.426155
[46,] -91.760744 -88.833624
[47,] -52.026022 -91.760744
[48,] -99.699162 -52.026022
[49,] 22.390478 -99.699162
[50,] -41.892266 22.390478
> z <- as.data.frame(dum1)
> z
lag(myerror, k = 1) myerror
1 60.521959 85.856374
2 -45.609982 60.521959
3 17.709353 -45.609982
4 -181.223659 17.709353
5 19.977176 -181.223659
6 -91.280613 19.977176
7 -45.689998 -91.280613
8 147.034075 -45.689998
9 -90.686608 147.034075
10 189.984027 -90.686608
11 120.250250 189.984027
12 -32.461252 120.250250
13 -78.174363 -32.461252
14 20.636723 -78.174363
15 25.624322 20.636723
16 -188.639454 25.624322
17 -66.449475 -188.639454
18 65.509498 -66.449475
19 16.312724 65.509498
20 -7.249921 16.312724
21 23.500947 -7.249921
22 43.516282 23.500947
23 -40.564182 43.516282
24 -1.142250 -40.564182
25 -5.450163 -1.142250
26 -17.887298 -5.450163
27 60.884764 -17.887298
28 -66.245814 60.884764
29 -156.351938 -66.245814
30 -43.692932 -156.351938
31 27.433545 -43.692932
32 -28.222446 27.433545
33 79.260493 -28.222446
34 29.127606 79.260493
35 -43.187717 29.127606
36 94.741286 -43.187717
37 -20.608911 94.741286
38 2.374403 -20.608911
39 109.977472 2.374403
40 150.368148 109.977472
41 25.784874 150.368148
42 130.164972 25.784874
43 38.662888 130.164972
44 17.426155 38.662888
45 -88.833624 17.426155
46 -91.760744 -88.833624
47 -52.026022 -91.760744
48 -99.699162 -52.026022
49 22.390478 -99.699162
50 -41.892266 22.390478
> 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/7h6hx1291646919.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/8xjn61291646919.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/9xjn61291646919.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/10xjn61291646919.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/116pwq1291646919.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/12zhdb1291646919.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/13n0sn1291646919.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/14909b1291646919.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/15ujph1291646919.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/16g1o41291646919.tab")
+ }
>
> try(system("convert tmp/1wo091291646919.ps tmp/1wo091291646919.png",intern=TRUE))
character(0)
> try(system("convert tmp/2wo091291646919.ps tmp/2wo091291646919.png",intern=TRUE))
character(0)
> try(system("convert tmp/3wo091291646919.ps tmp/3wo091291646919.png",intern=TRUE))
character(0)
> try(system("convert tmp/46xit1291646919.ps tmp/46xit1291646919.png",intern=TRUE))
character(0)
> try(system("convert tmp/56xit1291646919.ps tmp/56xit1291646919.png",intern=TRUE))
character(0)
> try(system("convert tmp/66xit1291646919.ps tmp/66xit1291646919.png",intern=TRUE))
character(0)
> try(system("convert tmp/7h6hx1291646919.ps tmp/7h6hx1291646919.png",intern=TRUE))
character(0)
> try(system("convert tmp/8xjn61291646919.ps tmp/8xjn61291646919.png",intern=TRUE))
character(0)
> try(system("convert tmp/9xjn61291646919.ps tmp/9xjn61291646919.png",intern=TRUE))
character(0)
> try(system("convert tmp/10xjn61291646919.ps tmp/10xjn61291646919.png",intern=TRUE))
character(0)
>
>
> proc.time()
user system elapsed
2.453 1.695 6.316