R version 2.12.1 (2010-12-16)
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(6.9
+ ,3.75
+ ,-16.16
+ ,-29.17
+ ,7.0
+ ,3.67
+ ,-14.17
+ ,-11.08
+ ,7.0
+ ,3.35
+ ,27.78
+ ,1.97
+ ,7.0
+ ,3.41
+ ,-15.92
+ ,-15.67
+ ,7.0
+ ,3.52
+ ,8.76
+ ,18.34
+ ,7.1
+ ,3.39
+ ,9.57
+ ,3.63
+ ,7.3
+ ,3.22
+ ,-1.64
+ ,68.84
+ ,7.6
+ ,3.1
+ ,15.75
+ ,-27.76
+ ,7.9
+ ,2.86
+ ,-8.31
+ ,-7.73
+ ,8.1
+ ,3.01
+ ,1.53
+ ,23.93
+ ,8.2
+ ,2.91
+ ,-1.34
+ ,-9.19
+ ,8.3
+ ,2.32
+ ,-39.82
+ ,0.91
+ ,8.5
+ ,2.57
+ ,-22.76
+ ,-27.81
+ ,8.5
+ ,2.46
+ ,16.42
+ ,8.83
+ ,8.5
+ ,2.27
+ ,7.50
+ ,-15.22
+ ,8.5
+ ,1.8
+ ,-5.17
+ ,-7.51
+ ,8.5
+ ,1.66
+ ,34.53
+ ,33.42
+ ,8.4
+ ,0.7
+ ,14.92
+ ,-5.72
+ ,8.3
+ ,0.62
+ ,1.45
+ ,89.58
+ ,8.2
+ ,0.26
+ ,2.05
+ ,-22.38
+ ,8.1
+ ,-0.12
+ ,-11.54
+ ,-16.70
+ ,8.0
+ ,-0.97
+ ,-6.86
+ ,12.22
+ ,8.1
+ ,-1.19
+ ,-1.09
+ ,18.06
+ ,8.1
+ ,-0.78
+ ,7.14
+ ,-10.81
+ ,8.0
+ ,-1.68
+ ,-10.35
+ ,-16.39
+ ,7.8
+ ,-1.1
+ ,17.18
+ ,10.93
+ ,7.7
+ ,-0.37
+ ,-6.67
+ ,-21.68
+ ,7.7
+ ,0.6
+ ,-0.13
+ ,-5.73
+ ,7.8
+ ,0.62
+ ,11.06
+ ,15.34
+ ,7.7
+ ,1.93
+ ,10.03
+ ,-5.33
+ ,7.5
+ ,2.32
+ ,-32.56
+ ,121.11
+ ,7.1
+ ,2.63
+ ,16.38
+ ,-26.09
+ ,7.0
+ ,3.14
+ ,2.19
+ ,-30.85
+ ,7.1
+ ,4.72
+ ,-9.34
+ ,9.74
+ ,7.3
+ ,5.46
+ ,16.01
+ ,15.35
+ ,7.4
+ ,5.39
+ ,-3.16
+ ,-11.96
+ ,7.3
+ ,5.91
+ ,-17.66
+ ,-24.71
+ ,6.9
+ ,5.8
+ ,16.98
+ ,3.91
+ ,6.7
+ ,5.21
+ ,-10.20
+ ,-19.13
+ ,6.7
+ ,4.15
+ ,3.60
+ ,6.38
+ ,6.8
+ ,4.39
+ ,-8.84
+ ,-0.31
+ ,6.9
+ ,3.64
+ ,2.17
+ ,-4.61
+ ,7.1
+ ,3.46
+ ,23.66
+ ,147.81
+ ,7.2
+ ,3.09
+ ,-9.86
+ ,-35.71
+ ,7.1
+ ,2.94
+ ,-24.89
+ ,-20.03
+ ,7.1
+ ,2.24
+ ,48.52
+ ,21.31
+ ,6.9
+ ,1.51
+ ,-16.80
+ ,9.29
+ ,7.0
+ ,1.12
+ ,6.47
+ ,-8.47
+ ,7.2
+ ,1.37
+ ,-15.74
+ ,-20.79
+ ,7.5
+ ,1.29
+ ,23.83
+ ,2.70
+ ,7.9
+ ,1.28
+ ,-2.00
+ ,-2.94
+ ,8.0
+ ,1.78
+ ,-17.81
+ ,-15.63
+ ,7.9
+ ,1.82
+ ,21.26
+ ,15.44
+ ,7.9
+ ,1.77
+ ,-10.92
+ ,-14.69
+ ,7.9
+ ,1.66
+ ,2.40
+ ,232.27
+ ,8.0
+ ,1.64
+ ,11.55
+ ,-46.11
+ ,8.0
+ ,1.49
+ ,-21.76
+ ,-17.21
+ ,8.0
+ ,1.21
+ ,5.75
+ ,15.33
+ ,8.0
+ ,1.22
+ ,5.92
+ ,5.06
+ ,7.9
+ ,1.63
+ ,2.42
+ ,-4.57
+ ,8.0
+ ,1.6
+ ,-17.19
+ ,-23.36
+ ,8.4
+ ,1.87
+ ,3.89
+ ,-11.39)
+ ,dim=c(4
+ ,62)
+ ,dimnames=list(c('Werkloosheid'
+ ,'inflatie'
+ ,'nieuwe_res_woningen'
+ ,'private_voertuigen')
+ ,1:62))
> y <- array(NA,dim=c(4,62),dimnames=list(c('Werkloosheid','inflatie','nieuwe_res_woningen','private_voertuigen'),1:62))
> 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'
> 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
Werkloosheid inflatie nieuwe_res_woningen private_voertuigen
1 6.9 3.75 -16.16 -29.17
2 7.0 3.67 -14.17 -11.08
3 7.0 3.35 27.78 1.97
4 7.0 3.41 -15.92 -15.67
5 7.0 3.52 8.76 18.34
6 7.1 3.39 9.57 3.63
7 7.3 3.22 -1.64 68.84
8 7.6 3.10 15.75 -27.76
9 7.9 2.86 -8.31 -7.73
10 8.1 3.01 1.53 23.93
11 8.2 2.91 -1.34 -9.19
12 8.3 2.32 -39.82 0.91
13 8.5 2.57 -22.76 -27.81
14 8.5 2.46 16.42 8.83
15 8.5 2.27 7.50 -15.22
16 8.5 1.80 -5.17 -7.51
17 8.5 1.66 34.53 33.42
18 8.4 0.70 14.92 -5.72
19 8.3 0.62 1.45 89.58
20 8.2 0.26 2.05 -22.38
21 8.1 -0.12 -11.54 -16.70
22 8.0 -0.97 -6.86 12.22
23 8.1 -1.19 -1.09 18.06
24 8.1 -0.78 7.14 -10.81
25 8.0 -1.68 -10.35 -16.39
26 7.8 -1.10 17.18 10.93
27 7.7 -0.37 -6.67 -21.68
28 7.7 0.60 -0.13 -5.73
29 7.8 0.62 11.06 15.34
30 7.7 1.93 10.03 -5.33
31 7.5 2.32 -32.56 121.11
32 7.1 2.63 16.38 -26.09
33 7.0 3.14 2.19 -30.85
34 7.1 4.72 -9.34 9.74
35 7.3 5.46 16.01 15.35
36 7.4 5.39 -3.16 -11.96
37 7.3 5.91 -17.66 -24.71
38 6.9 5.80 16.98 3.91
39 6.7 5.21 -10.20 -19.13
40 6.7 4.15 3.60 6.38
41 6.8 4.39 -8.84 -0.31
42 6.9 3.64 2.17 -4.61
43 7.1 3.46 23.66 147.81
44 7.2 3.09 -9.86 -35.71
45 7.1 2.94 -24.89 -20.03
46 7.1 2.24 48.52 21.31
47 6.9 1.51 -16.80 9.29
48 7.0 1.12 6.47 -8.47
49 7.2 1.37 -15.74 -20.79
50 7.5 1.29 23.83 2.70
51 7.9 1.28 -2.00 -2.94
52 8.0 1.78 -17.81 -15.63
53 7.9 1.82 21.26 15.44
54 7.9 1.77 -10.92 -14.69
55 7.9 1.66 2.40 232.27
56 8.0 1.64 11.55 -46.11
57 8.0 1.49 -21.76 -17.21
58 8.0 1.21 5.75 15.33
59 8.0 1.22 5.92 5.06
60 7.9 1.63 2.42 -4.57
61 8.0 1.60 -17.19 -23.36
62 8.4 1.87 3.89 -11.39
> k <- length(x[1,])
> df <- as.data.frame(x)
> (mylm <- lm(df))
Call:
lm(formula = df)
Coefficients:
(Intercept) inflatie nieuwe_res_woningen
8.0399545 -0.1831576 -0.0010197
private_voertuigen
0.0001561
> (mysum <- summary(mylm))
Call:
lm(formula = df)
Residuals:
Min 1Q Median 3Q Max
-0.88197 -0.39213 -0.06739 0.27310 0.92598
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 8.0399545 0.0954661 84.218 < 2e-16 ***
inflatie -0.1831576 0.0340699 -5.376 1.42e-06 ***
nieuwe_res_woningen -0.0010197 0.0036941 -0.276 0.784
private_voertuigen 0.0001561 0.0013474 0.116 0.908
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 0.4601 on 58 degrees of freedom
Multiple R-squared: 0.333, Adjusted R-squared: 0.2985
F-statistic: 9.652 on 3 and 58 DF, p-value: 2.912e-05
> 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.001843742 0.003687484 0.998156258
[2,] 0.017690237 0.035380474 0.982309763
[3,] 0.005418590 0.010837181 0.994581410
[4,] 0.028550496 0.057100991 0.971449504
[5,] 0.022484814 0.044969629 0.977515186
[6,] 0.045845912 0.091691824 0.954154088
[7,] 0.039536161 0.079072323 0.960463839
[8,] 0.028882697 0.057765394 0.971117303
[9,] 0.028928176 0.057856352 0.971071824
[10,] 0.123375205 0.246750411 0.876624795
[11,] 0.186017924 0.372035848 0.813982076
[12,] 0.709121472 0.581757056 0.290878528
[13,] 0.790770852 0.418458297 0.209229148
[14,] 0.914420199 0.171159602 0.085579801
[15,] 0.948765097 0.102469806 0.051234903
[16,] 0.968286800 0.063426401 0.031713200
[17,] 0.966083041 0.067833918 0.033916959
[18,] 0.956256839 0.087486322 0.043743161
[19,] 0.951087266 0.097825468 0.048912734
[20,] 0.947777850 0.104444301 0.052222150
[21,] 0.940346024 0.119307952 0.059653976
[22,] 0.920365040 0.159269920 0.079634960
[23,] 0.890539859 0.218920282 0.109460141
[24,] 0.852662136 0.294675728 0.147337864
[25,] 0.813268472 0.373463056 0.186731528
[26,] 0.811879804 0.376240392 0.188120196
[27,] 0.810917367 0.378165266 0.189082633
[28,] 0.758477447 0.483045106 0.241522553
[29,] 0.725290017 0.549419967 0.274709983
[30,] 0.718162308 0.563675384 0.281837692
[31,] 0.747860667 0.504278667 0.252139333
[32,] 0.738188442 0.523623117 0.261811558
[33,] 0.704235394 0.591529213 0.295764606
[34,] 0.683806161 0.632387679 0.316193839
[35,] 0.632005951 0.735988099 0.367994049
[36,] 0.585473727 0.829052546 0.414526273
[37,] 0.515744740 0.968510520 0.484255260
[38,] 0.445463664 0.890927328 0.554536336
[39,] 0.510659949 0.978680102 0.489340051
[40,] 0.742432965 0.515134069 0.257567035
[41,] 0.929882393 0.140235214 0.070117607
[42,] 0.972264290 0.055471421 0.027735710
[43,] 0.997967391 0.004065218 0.002032609
[44,] 0.999240639 0.001518722 0.000759361
[45,] 0.997578424 0.004843152 0.002421576
[46,] 0.992522665 0.014954671 0.007477335
[47,] 0.986363310 0.027273380 0.013636690
[48,] 0.972903524 0.054192952 0.027096476
[49,] 0.945497946 0.109004107 0.054502054
> postscript(file="/var/www/rcomp/tmp/1795x1321903651.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/2ddd31321903651.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/3c1771321903651.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/4i1m21321903651.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/539av1321903651.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 = 62
Frequency = 1
1 2 3 4 5 6
-0.46503770 -0.38048508 -0.39835784 -0.42917395 -0.38917042 -0.30985871
7 8 9 10 11 12
-0.16260539 0.14822717 0.37660954 0.60917443 0.69310237 0.64422617
13 14 15 16 17 18
0.91189427 0.92597768 0.88583663 0.78562989 0.79407910 0.50436211
19 20 21 22 23 24
0.36109805 0.21325040 0.02890665 -0.22651979 -0.16184265 -0.07384953
25 26 27 28 29 30
-0.35565420 -0.42561624 -0.41113961 -0.22929799 -0.11751392 0.02459893
31 32 33 34 35 36
-0.16713467 -0.43747520 -0.45779078 -0.08649470 0.27401461 0.34590983
37 38 39 40 41 42
0.32835700 -0.06093692 -0.39311770 -0.57717560 -0.44485804 -0.47032852
43 44 45 46 47 48
-0.30517757 -0.27847693 -0.42372378 -0.48353402 -0.88196701 -0.82689856
49 50 51 52 53 54
-0.60183266 -0.27980411 0.09290687 0.27034577 0.21266014 0.17539292
55 56 57 58 59 60
0.13027629 0.27939899 0.21344904 0.18513622 0.18874432 0.16177339
61 62
0.23921627 0.70829473
> postscript(file="/var/www/rcomp/tmp/6rpcz1321903651.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 = 62
Frequency = 1
lag(myerror, k = 1) myerror
0 -0.46503770 NA
1 -0.38048508 -0.46503770
2 -0.39835784 -0.38048508
3 -0.42917395 -0.39835784
4 -0.38917042 -0.42917395
5 -0.30985871 -0.38917042
6 -0.16260539 -0.30985871
7 0.14822717 -0.16260539
8 0.37660954 0.14822717
9 0.60917443 0.37660954
10 0.69310237 0.60917443
11 0.64422617 0.69310237
12 0.91189427 0.64422617
13 0.92597768 0.91189427
14 0.88583663 0.92597768
15 0.78562989 0.88583663
16 0.79407910 0.78562989
17 0.50436211 0.79407910
18 0.36109805 0.50436211
19 0.21325040 0.36109805
20 0.02890665 0.21325040
21 -0.22651979 0.02890665
22 -0.16184265 -0.22651979
23 -0.07384953 -0.16184265
24 -0.35565420 -0.07384953
25 -0.42561624 -0.35565420
26 -0.41113961 -0.42561624
27 -0.22929799 -0.41113961
28 -0.11751392 -0.22929799
29 0.02459893 -0.11751392
30 -0.16713467 0.02459893
31 -0.43747520 -0.16713467
32 -0.45779078 -0.43747520
33 -0.08649470 -0.45779078
34 0.27401461 -0.08649470
35 0.34590983 0.27401461
36 0.32835700 0.34590983
37 -0.06093692 0.32835700
38 -0.39311770 -0.06093692
39 -0.57717560 -0.39311770
40 -0.44485804 -0.57717560
41 -0.47032852 -0.44485804
42 -0.30517757 -0.47032852
43 -0.27847693 -0.30517757
44 -0.42372378 -0.27847693
45 -0.48353402 -0.42372378
46 -0.88196701 -0.48353402
47 -0.82689856 -0.88196701
48 -0.60183266 -0.82689856
49 -0.27980411 -0.60183266
50 0.09290687 -0.27980411
51 0.27034577 0.09290687
52 0.21266014 0.27034577
53 0.17539292 0.21266014
54 0.13027629 0.17539292
55 0.27939899 0.13027629
56 0.21344904 0.27939899
57 0.18513622 0.21344904
58 0.18874432 0.18513622
59 0.16177339 0.18874432
60 0.23921627 0.16177339
61 0.70829473 0.23921627
62 NA 0.70829473
> dum1 <- dum[2:length(myerror),]
> dum1
lag(myerror, k = 1) myerror
[1,] -0.38048508 -0.46503770
[2,] -0.39835784 -0.38048508
[3,] -0.42917395 -0.39835784
[4,] -0.38917042 -0.42917395
[5,] -0.30985871 -0.38917042
[6,] -0.16260539 -0.30985871
[7,] 0.14822717 -0.16260539
[8,] 0.37660954 0.14822717
[9,] 0.60917443 0.37660954
[10,] 0.69310237 0.60917443
[11,] 0.64422617 0.69310237
[12,] 0.91189427 0.64422617
[13,] 0.92597768 0.91189427
[14,] 0.88583663 0.92597768
[15,] 0.78562989 0.88583663
[16,] 0.79407910 0.78562989
[17,] 0.50436211 0.79407910
[18,] 0.36109805 0.50436211
[19,] 0.21325040 0.36109805
[20,] 0.02890665 0.21325040
[21,] -0.22651979 0.02890665
[22,] -0.16184265 -0.22651979
[23,] -0.07384953 -0.16184265
[24,] -0.35565420 -0.07384953
[25,] -0.42561624 -0.35565420
[26,] -0.41113961 -0.42561624
[27,] -0.22929799 -0.41113961
[28,] -0.11751392 -0.22929799
[29,] 0.02459893 -0.11751392
[30,] -0.16713467 0.02459893
[31,] -0.43747520 -0.16713467
[32,] -0.45779078 -0.43747520
[33,] -0.08649470 -0.45779078
[34,] 0.27401461 -0.08649470
[35,] 0.34590983 0.27401461
[36,] 0.32835700 0.34590983
[37,] -0.06093692 0.32835700
[38,] -0.39311770 -0.06093692
[39,] -0.57717560 -0.39311770
[40,] -0.44485804 -0.57717560
[41,] -0.47032852 -0.44485804
[42,] -0.30517757 -0.47032852
[43,] -0.27847693 -0.30517757
[44,] -0.42372378 -0.27847693
[45,] -0.48353402 -0.42372378
[46,] -0.88196701 -0.48353402
[47,] -0.82689856 -0.88196701
[48,] -0.60183266 -0.82689856
[49,] -0.27980411 -0.60183266
[50,] 0.09290687 -0.27980411
[51,] 0.27034577 0.09290687
[52,] 0.21266014 0.27034577
[53,] 0.17539292 0.21266014
[54,] 0.13027629 0.17539292
[55,] 0.27939899 0.13027629
[56,] 0.21344904 0.27939899
[57,] 0.18513622 0.21344904
[58,] 0.18874432 0.18513622
[59,] 0.16177339 0.18874432
[60,] 0.23921627 0.16177339
[61,] 0.70829473 0.23921627
> z <- as.data.frame(dum1)
> z
lag(myerror, k = 1) myerror
1 -0.38048508 -0.46503770
2 -0.39835784 -0.38048508
3 -0.42917395 -0.39835784
4 -0.38917042 -0.42917395
5 -0.30985871 -0.38917042
6 -0.16260539 -0.30985871
7 0.14822717 -0.16260539
8 0.37660954 0.14822717
9 0.60917443 0.37660954
10 0.69310237 0.60917443
11 0.64422617 0.69310237
12 0.91189427 0.64422617
13 0.92597768 0.91189427
14 0.88583663 0.92597768
15 0.78562989 0.88583663
16 0.79407910 0.78562989
17 0.50436211 0.79407910
18 0.36109805 0.50436211
19 0.21325040 0.36109805
20 0.02890665 0.21325040
21 -0.22651979 0.02890665
22 -0.16184265 -0.22651979
23 -0.07384953 -0.16184265
24 -0.35565420 -0.07384953
25 -0.42561624 -0.35565420
26 -0.41113961 -0.42561624
27 -0.22929799 -0.41113961
28 -0.11751392 -0.22929799
29 0.02459893 -0.11751392
30 -0.16713467 0.02459893
31 -0.43747520 -0.16713467
32 -0.45779078 -0.43747520
33 -0.08649470 -0.45779078
34 0.27401461 -0.08649470
35 0.34590983 0.27401461
36 0.32835700 0.34590983
37 -0.06093692 0.32835700
38 -0.39311770 -0.06093692
39 -0.57717560 -0.39311770
40 -0.44485804 -0.57717560
41 -0.47032852 -0.44485804
42 -0.30517757 -0.47032852
43 -0.27847693 -0.30517757
44 -0.42372378 -0.27847693
45 -0.48353402 -0.42372378
46 -0.88196701 -0.48353402
47 -0.82689856 -0.88196701
48 -0.60183266 -0.82689856
49 -0.27980411 -0.60183266
50 0.09290687 -0.27980411
51 0.27034577 0.09290687
52 0.21266014 0.27034577
53 0.17539292 0.21266014
54 0.13027629 0.17539292
55 0.27939899 0.13027629
56 0.21344904 0.27939899
57 0.18513622 0.21344904
58 0.18874432 0.18513622
59 0.16177339 0.18874432
60 0.23921627 0.16177339
61 0.70829473 0.23921627
> 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/706ui1321903651.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/89fag1321903651.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/9pge31321903651.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/10ts661321903651.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/11teec1321903651.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/12f0ej1321903652.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/13aof11321903652.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/14er1e1321903652.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/15yww71321903652.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/16b7o11321903652.tab")
+ }
>
> try(system("convert tmp/1795x1321903651.ps tmp/1795x1321903651.png",intern=TRUE))
character(0)
> try(system("convert tmp/2ddd31321903651.ps tmp/2ddd31321903651.png",intern=TRUE))
character(0)
> try(system("convert tmp/3c1771321903651.ps tmp/3c1771321903651.png",intern=TRUE))
character(0)
> try(system("convert tmp/4i1m21321903651.ps tmp/4i1m21321903651.png",intern=TRUE))
character(0)
> try(system("convert tmp/539av1321903651.ps tmp/539av1321903651.png",intern=TRUE))
character(0)
> try(system("convert tmp/6rpcz1321903651.ps tmp/6rpcz1321903651.png",intern=TRUE))
character(0)
> try(system("convert tmp/706ui1321903651.ps tmp/706ui1321903651.png",intern=TRUE))
character(0)
> try(system("convert tmp/89fag1321903651.ps tmp/89fag1321903651.png",intern=TRUE))
character(0)
> try(system("convert tmp/9pge31321903651.ps tmp/9pge31321903651.png",intern=TRUE))
character(0)
> try(system("convert tmp/10ts661321903651.ps tmp/10ts661321903651.png",intern=TRUE))
character(0)
>
>
> proc.time()
user system elapsed
4.228 0.660 4.873