R version 2.13.0 (2011-04-13)
Copyright (C) 2011 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(19.785,18.479,10.698,31.956,29.506,34.506,27.165,26.736,23.691,18.157,17.328,18.205,20.995,17.382,9.367,31.124,26.551,30.651,25.859,25.100,25.778,20.418,18.688,20.424,24.776,19.814,12.738,31.566,30.111,30.019,31.934,25.826,26.835,20.205,17.789,20.520,22.518,15.572,11.509,25.447,24.090,27.786,26.195,20.516,22.759,19.028,16.971,20.036,22.485,18.730,14.538,27.561,25.985,34.670,32.066,27.186,29.586,21.359,21.553,19.573,24.256,22.380,16.167,27.297,28.287,33.474,28.229,28.785,25.597,18.130,20.198,22.849,23.118),dim=c(1,73),dimnames=list(c('Inschrijvingen'),1:73))
> y <- array(NA,dim=c(1,73),dimnames=list(c('Inschrijvingen'),1:73))
> 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 = 'Include Monthly 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
Inschrijvingen M1 M2 M3 M4 M5 M6 M7 M8 M9 M10 M11 t
1 19.785 1 0 0 0 0 0 0 0 0 0 0 1
2 18.479 0 1 0 0 0 0 0 0 0 0 0 2
3 10.698 0 0 1 0 0 0 0 0 0 0 0 3
4 31.956 0 0 0 1 0 0 0 0 0 0 0 4
5 29.506 0 0 0 0 1 0 0 0 0 0 0 5
6 34.506 0 0 0 0 0 1 0 0 0 0 0 6
7 27.165 0 0 0 0 0 0 1 0 0 0 0 7
8 26.736 0 0 0 0 0 0 0 1 0 0 0 8
9 23.691 0 0 0 0 0 0 0 0 1 0 0 9
10 18.157 0 0 0 0 0 0 0 0 0 1 0 10
11 17.328 0 0 0 0 0 0 0 0 0 0 1 11
12 18.205 0 0 0 0 0 0 0 0 0 0 0 12
13 20.995 1 0 0 0 0 0 0 0 0 0 0 13
14 17.382 0 1 0 0 0 0 0 0 0 0 0 14
15 9.367 0 0 1 0 0 0 0 0 0 0 0 15
16 31.124 0 0 0 1 0 0 0 0 0 0 0 16
17 26.551 0 0 0 0 1 0 0 0 0 0 0 17
18 30.651 0 0 0 0 0 1 0 0 0 0 0 18
19 25.859 0 0 0 0 0 0 1 0 0 0 0 19
20 25.100 0 0 0 0 0 0 0 1 0 0 0 20
21 25.778 0 0 0 0 0 0 0 0 1 0 0 21
22 20.418 0 0 0 0 0 0 0 0 0 1 0 22
23 18.688 0 0 0 0 0 0 0 0 0 0 1 23
24 20.424 0 0 0 0 0 0 0 0 0 0 0 24
25 24.776 1 0 0 0 0 0 0 0 0 0 0 25
26 19.814 0 1 0 0 0 0 0 0 0 0 0 26
27 12.738 0 0 1 0 0 0 0 0 0 0 0 27
28 31.566 0 0 0 1 0 0 0 0 0 0 0 28
29 30.111 0 0 0 0 1 0 0 0 0 0 0 29
30 30.019 0 0 0 0 0 1 0 0 0 0 0 30
31 31.934 0 0 0 0 0 0 1 0 0 0 0 31
32 25.826 0 0 0 0 0 0 0 1 0 0 0 32
33 26.835 0 0 0 0 0 0 0 0 1 0 0 33
34 20.205 0 0 0 0 0 0 0 0 0 1 0 34
35 17.789 0 0 0 0 0 0 0 0 0 0 1 35
36 20.520 0 0 0 0 0 0 0 0 0 0 0 36
37 22.518 1 0 0 0 0 0 0 0 0 0 0 37
38 15.572 0 1 0 0 0 0 0 0 0 0 0 38
39 11.509 0 0 1 0 0 0 0 0 0 0 0 39
40 25.447 0 0 0 1 0 0 0 0 0 0 0 40
41 24.090 0 0 0 0 1 0 0 0 0 0 0 41
42 27.786 0 0 0 0 0 1 0 0 0 0 0 42
43 26.195 0 0 0 0 0 0 1 0 0 0 0 43
44 20.516 0 0 0 0 0 0 0 1 0 0 0 44
45 22.759 0 0 0 0 0 0 0 0 1 0 0 45
46 19.028 0 0 0 0 0 0 0 0 0 1 0 46
47 16.971 0 0 0 0 0 0 0 0 0 0 1 47
48 20.036 0 0 0 0 0 0 0 0 0 0 0 48
49 22.485 1 0 0 0 0 0 0 0 0 0 0 49
50 18.730 0 1 0 0 0 0 0 0 0 0 0 50
51 14.538 0 0 1 0 0 0 0 0 0 0 0 51
52 27.561 0 0 0 1 0 0 0 0 0 0 0 52
53 25.985 0 0 0 0 1 0 0 0 0 0 0 53
54 34.670 0 0 0 0 0 1 0 0 0 0 0 54
55 32.066 0 0 0 0 0 0 1 0 0 0 0 55
56 27.186 0 0 0 0 0 0 0 1 0 0 0 56
57 29.586 0 0 0 0 0 0 0 0 1 0 0 57
58 21.359 0 0 0 0 0 0 0 0 0 1 0 58
59 21.553 0 0 0 0 0 0 0 0 0 0 1 59
60 19.573 0 0 0 0 0 0 0 0 0 0 0 60
61 24.256 1 0 0 0 0 0 0 0 0 0 0 61
62 22.380 0 1 0 0 0 0 0 0 0 0 0 62
63 16.167 0 0 1 0 0 0 0 0 0 0 0 63
64 27.297 0 0 0 1 0 0 0 0 0 0 0 64
65 28.287 0 0 0 0 1 0 0 0 0 0 0 65
66 33.474 0 0 0 0 0 1 0 0 0 0 0 66
67 28.229 0 0 0 0 0 0 1 0 0 0 0 67
68 28.785 0 0 0 0 0 0 0 1 0 0 0 68
69 25.597 0 0 0 0 0 0 0 0 1 0 0 69
70 18.130 0 0 0 0 0 0 0 0 0 1 0 70
71 20.198 0 0 0 0 0 0 0 0 0 0 1 71
72 22.849 0 0 0 0 0 0 0 0 0 0 0 72
73 23.118 1 0 0 0 0 0 0 0 0 0 0 73
> k <- length(x[1,])
> df <- as.data.frame(x)
> (mylm <- lm(df))
Call:
lm(formula = df)
Coefficients:
(Intercept) M1 M2 M3 M4 M5
19.23421 2.41707 -1.29557 -7.54351 9.08755 7.32610
M6 M7 M8 M9 M10 M11
11.73083 8.42988 5.52211 5.51366 -0.66911 -1.48872
t
0.02461
> (mysum <- summary(mylm))
Call:
lm(formula = df)
Residuals:
Min 1Q Median 3Q Max
-5.3232 -1.1415 -0.1485 1.3665 3.5358
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 19.23421 1.06041 18.139 < 2e-16 ***
M1 2.41707 1.25308 1.929 0.0585 .
M2 -1.29557 1.30485 -0.993 0.3248
M3 -7.54351 1.30369 -5.786 2.79e-07 ***
M4 9.08755 1.30265 6.976 2.75e-09 ***
M5 7.32610 1.30173 5.628 5.09e-07 ***
M6 11.73083 1.30094 9.017 9.23e-13 ***
M7 8.42988 1.30026 6.483 1.89e-08 ***
M8 5.52211 1.29971 4.249 7.61e-05 ***
M9 5.51366 1.29928 4.244 7.74e-05 ***
M10 -0.66911 1.29897 -0.515 0.6084
M11 -1.48872 1.29879 -1.146 0.2562
t 0.02461 0.01262 1.949 0.0559 .
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 2.249 on 60 degrees of freedom
Multiple R-squared: 0.8734, Adjusted R-squared: 0.848
F-statistic: 34.48 on 12 and 60 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.06050013 0.12100026 0.9394999
[2,] 0.06997417 0.13994835 0.9300258
[3,] 0.08614053 0.17228106 0.9138595
[4,] 0.04203067 0.08406134 0.9579693
[5,] 0.01796118 0.03592236 0.9820388
[6,] 0.03963538 0.07927077 0.9603646
[7,] 0.05474203 0.10948407 0.9452580
[8,] 0.04136374 0.08272748 0.9586363
[9,] 0.03933302 0.07866604 0.9606670
[10,] 0.10894181 0.21788361 0.8910582
[11,] 0.07980013 0.15960026 0.9201999
[12,] 0.06136092 0.12272183 0.9386391
[13,] 0.06984477 0.13968953 0.9301552
[14,] 0.08320654 0.16641307 0.9167935
[15,] 0.10507871 0.21015743 0.8949213
[16,] 0.25978499 0.51956999 0.7402150
[17,] 0.22617926 0.45235852 0.7738207
[18,] 0.21251664 0.42503327 0.7874834
[19,] 0.19887584 0.39775168 0.8011242
[20,] 0.15781869 0.31563739 0.8421813
[21,] 0.13710763 0.27421526 0.8628924
[22,] 0.11861836 0.23723672 0.8813816
[23,] 0.18486514 0.36973027 0.8151349
[24,] 0.14151774 0.28303549 0.8584823
[25,] 0.32946708 0.65893415 0.6705329
[26,] 0.37493128 0.74986256 0.6250687
[27,] 0.47319757 0.94639515 0.5268024
[28,] 0.43114856 0.86229711 0.5688514
[29,] 0.73007062 0.53985876 0.2699294
[30,] 0.79186114 0.41627772 0.2081389
[31,] 0.71916327 0.56167347 0.2808367
[32,] 0.75718357 0.48563285 0.2428164
[33,] 0.69262416 0.61475169 0.3073758
[34,] 0.63045928 0.73908145 0.3695407
[35,] 0.70332193 0.59335613 0.2966781
[36,] 0.71166888 0.57666223 0.2883311
[37,] 0.62061173 0.75877654 0.3793883
[38,] 0.64086706 0.71826589 0.3591329
[39,] 0.56815297 0.86369406 0.4318470
[40,] 0.58201152 0.83597696 0.4179885
[41,] 0.54498196 0.91003608 0.4550180
[42,] 0.55640028 0.88719944 0.4435997
> postscript(file="/var/wessaorg/rcomp/tmp/16qez1322582577.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/wessaorg/rcomp/tmp/28gxb1322582577.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/wessaorg/rcomp/tmp/33vhf1322582577.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/wessaorg/rcomp/tmp/4i1eq1322582577.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/wessaorg/rcomp/tmp/5zt741322582578.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 = 73
Frequency = 1
1 2 3 4 5 6
-1.89089116 0.49113832 -1.06652834 3.53580499 2.82263832 3.39330499
7 8 9 10 11 12
-0.67136168 1.78280499 -1.27836168 -0.65419501 -0.68819501 -1.32452834
13 14 15 16 17 18
-0.97621315 -0.90118367 -2.69285034 2.40848299 -0.42768367 -0.75701701
19 20 21 22 23 24
-2.27268367 -0.14851701 0.51331633 1.31148299 0.37648299 0.59914966
25 26 27 28 29 30
2.50946485 1.23549433 0.38282766 2.55516100 2.83699433 -1.68433900
31 32 33 34 35 36
3.50699433 0.28216100 1.27499433 0.80316100 -0.81783900 0.39982766
37 38 39 40 41 42
-0.04385714 -3.30182766 -1.14149433 -3.85916100 -3.47932766 -4.21266100
43 44 45 46 47 48
-2.52732766 -5.32316100 -3.09632766 -0.66916100 -1.93116100 -0.37949433
49 50 51 52 53 54
-0.37217914 -0.43914966 1.59218367 -2.04048299 -1.87964966 2.37601701
55 56 57 58 59 60
3.04835034 1.05151701 3.43535034 1.36651701 2.35551701 -1.13781633
61 62 63 64 65 66
1.10349887 2.91552834 2.92586168 -2.59980499 0.12702834 0.88469501
67 68 69 70 71 72
-1.08397166 2.35519501 -0.84897166 -2.15780499 0.70519501 1.84286168
73
-0.32982313
> postscript(file="/var/wessaorg/rcomp/tmp/6s7ii1322582578.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 = 73
Frequency = 1
lag(myerror, k = 1) myerror
0 -1.89089116 NA
1 0.49113832 -1.89089116
2 -1.06652834 0.49113832
3 3.53580499 -1.06652834
4 2.82263832 3.53580499
5 3.39330499 2.82263832
6 -0.67136168 3.39330499
7 1.78280499 -0.67136168
8 -1.27836168 1.78280499
9 -0.65419501 -1.27836168
10 -0.68819501 -0.65419501
11 -1.32452834 -0.68819501
12 -0.97621315 -1.32452834
13 -0.90118367 -0.97621315
14 -2.69285034 -0.90118367
15 2.40848299 -2.69285034
16 -0.42768367 2.40848299
17 -0.75701701 -0.42768367
18 -2.27268367 -0.75701701
19 -0.14851701 -2.27268367
20 0.51331633 -0.14851701
21 1.31148299 0.51331633
22 0.37648299 1.31148299
23 0.59914966 0.37648299
24 2.50946485 0.59914966
25 1.23549433 2.50946485
26 0.38282766 1.23549433
27 2.55516100 0.38282766
28 2.83699433 2.55516100
29 -1.68433900 2.83699433
30 3.50699433 -1.68433900
31 0.28216100 3.50699433
32 1.27499433 0.28216100
33 0.80316100 1.27499433
34 -0.81783900 0.80316100
35 0.39982766 -0.81783900
36 -0.04385714 0.39982766
37 -3.30182766 -0.04385714
38 -1.14149433 -3.30182766
39 -3.85916100 -1.14149433
40 -3.47932766 -3.85916100
41 -4.21266100 -3.47932766
42 -2.52732766 -4.21266100
43 -5.32316100 -2.52732766
44 -3.09632766 -5.32316100
45 -0.66916100 -3.09632766
46 -1.93116100 -0.66916100
47 -0.37949433 -1.93116100
48 -0.37217914 -0.37949433
49 -0.43914966 -0.37217914
50 1.59218367 -0.43914966
51 -2.04048299 1.59218367
52 -1.87964966 -2.04048299
53 2.37601701 -1.87964966
54 3.04835034 2.37601701
55 1.05151701 3.04835034
56 3.43535034 1.05151701
57 1.36651701 3.43535034
58 2.35551701 1.36651701
59 -1.13781633 2.35551701
60 1.10349887 -1.13781633
61 2.91552834 1.10349887
62 2.92586168 2.91552834
63 -2.59980499 2.92586168
64 0.12702834 -2.59980499
65 0.88469501 0.12702834
66 -1.08397166 0.88469501
67 2.35519501 -1.08397166
68 -0.84897166 2.35519501
69 -2.15780499 -0.84897166
70 0.70519501 -2.15780499
71 1.84286168 0.70519501
72 -0.32982313 1.84286168
73 NA -0.32982313
> dum1 <- dum[2:length(myerror),]
> dum1
lag(myerror, k = 1) myerror
[1,] 0.49113832 -1.89089116
[2,] -1.06652834 0.49113832
[3,] 3.53580499 -1.06652834
[4,] 2.82263832 3.53580499
[5,] 3.39330499 2.82263832
[6,] -0.67136168 3.39330499
[7,] 1.78280499 -0.67136168
[8,] -1.27836168 1.78280499
[9,] -0.65419501 -1.27836168
[10,] -0.68819501 -0.65419501
[11,] -1.32452834 -0.68819501
[12,] -0.97621315 -1.32452834
[13,] -0.90118367 -0.97621315
[14,] -2.69285034 -0.90118367
[15,] 2.40848299 -2.69285034
[16,] -0.42768367 2.40848299
[17,] -0.75701701 -0.42768367
[18,] -2.27268367 -0.75701701
[19,] -0.14851701 -2.27268367
[20,] 0.51331633 -0.14851701
[21,] 1.31148299 0.51331633
[22,] 0.37648299 1.31148299
[23,] 0.59914966 0.37648299
[24,] 2.50946485 0.59914966
[25,] 1.23549433 2.50946485
[26,] 0.38282766 1.23549433
[27,] 2.55516100 0.38282766
[28,] 2.83699433 2.55516100
[29,] -1.68433900 2.83699433
[30,] 3.50699433 -1.68433900
[31,] 0.28216100 3.50699433
[32,] 1.27499433 0.28216100
[33,] 0.80316100 1.27499433
[34,] -0.81783900 0.80316100
[35,] 0.39982766 -0.81783900
[36,] -0.04385714 0.39982766
[37,] -3.30182766 -0.04385714
[38,] -1.14149433 -3.30182766
[39,] -3.85916100 -1.14149433
[40,] -3.47932766 -3.85916100
[41,] -4.21266100 -3.47932766
[42,] -2.52732766 -4.21266100
[43,] -5.32316100 -2.52732766
[44,] -3.09632766 -5.32316100
[45,] -0.66916100 -3.09632766
[46,] -1.93116100 -0.66916100
[47,] -0.37949433 -1.93116100
[48,] -0.37217914 -0.37949433
[49,] -0.43914966 -0.37217914
[50,] 1.59218367 -0.43914966
[51,] -2.04048299 1.59218367
[52,] -1.87964966 -2.04048299
[53,] 2.37601701 -1.87964966
[54,] 3.04835034 2.37601701
[55,] 1.05151701 3.04835034
[56,] 3.43535034 1.05151701
[57,] 1.36651701 3.43535034
[58,] 2.35551701 1.36651701
[59,] -1.13781633 2.35551701
[60,] 1.10349887 -1.13781633
[61,] 2.91552834 1.10349887
[62,] 2.92586168 2.91552834
[63,] -2.59980499 2.92586168
[64,] 0.12702834 -2.59980499
[65,] 0.88469501 0.12702834
[66,] -1.08397166 0.88469501
[67,] 2.35519501 -1.08397166
[68,] -0.84897166 2.35519501
[69,] -2.15780499 -0.84897166
[70,] 0.70519501 -2.15780499
[71,] 1.84286168 0.70519501
[72,] -0.32982313 1.84286168
> z <- as.data.frame(dum1)
> z
lag(myerror, k = 1) myerror
1 0.49113832 -1.89089116
2 -1.06652834 0.49113832
3 3.53580499 -1.06652834
4 2.82263832 3.53580499
5 3.39330499 2.82263832
6 -0.67136168 3.39330499
7 1.78280499 -0.67136168
8 -1.27836168 1.78280499
9 -0.65419501 -1.27836168
10 -0.68819501 -0.65419501
11 -1.32452834 -0.68819501
12 -0.97621315 -1.32452834
13 -0.90118367 -0.97621315
14 -2.69285034 -0.90118367
15 2.40848299 -2.69285034
16 -0.42768367 2.40848299
17 -0.75701701 -0.42768367
18 -2.27268367 -0.75701701
19 -0.14851701 -2.27268367
20 0.51331633 -0.14851701
21 1.31148299 0.51331633
22 0.37648299 1.31148299
23 0.59914966 0.37648299
24 2.50946485 0.59914966
25 1.23549433 2.50946485
26 0.38282766 1.23549433
27 2.55516100 0.38282766
28 2.83699433 2.55516100
29 -1.68433900 2.83699433
30 3.50699433 -1.68433900
31 0.28216100 3.50699433
32 1.27499433 0.28216100
33 0.80316100 1.27499433
34 -0.81783900 0.80316100
35 0.39982766 -0.81783900
36 -0.04385714 0.39982766
37 -3.30182766 -0.04385714
38 -1.14149433 -3.30182766
39 -3.85916100 -1.14149433
40 -3.47932766 -3.85916100
41 -4.21266100 -3.47932766
42 -2.52732766 -4.21266100
43 -5.32316100 -2.52732766
44 -3.09632766 -5.32316100
45 -0.66916100 -3.09632766
46 -1.93116100 -0.66916100
47 -0.37949433 -1.93116100
48 -0.37217914 -0.37949433
49 -0.43914966 -0.37217914
50 1.59218367 -0.43914966
51 -2.04048299 1.59218367
52 -1.87964966 -2.04048299
53 2.37601701 -1.87964966
54 3.04835034 2.37601701
55 1.05151701 3.04835034
56 3.43535034 1.05151701
57 1.36651701 3.43535034
58 2.35551701 1.36651701
59 -1.13781633 2.35551701
60 1.10349887 -1.13781633
61 2.91552834 1.10349887
62 2.92586168 2.91552834
63 -2.59980499 2.92586168
64 0.12702834 -2.59980499
65 0.88469501 0.12702834
66 -1.08397166 0.88469501
67 2.35519501 -1.08397166
68 -0.84897166 2.35519501
69 -2.15780499 -0.84897166
70 0.70519501 -2.15780499
71 1.84286168 0.70519501
72 -0.32982313 1.84286168
> 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/wessaorg/rcomp/tmp/7pmn91322582578.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/wessaorg/rcomp/tmp/8vlm21322582578.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/wessaorg/rcomp/tmp/9b01k1322582578.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/wessaorg/rcomp/tmp/10hyvk1322582578.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/wessaorg/rcomp/createtable file can be downloaded at http://www.wessa.net/cretab
> load(file="/var/wessaorg/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/wessaorg/rcomp/tmp/11udzc1322582578.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/wessaorg/rcomp/tmp/12fd2f1322582578.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/wessaorg/rcomp/tmp/1379jz1322582578.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/wessaorg/rcomp/tmp/14h2h61322582578.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/wessaorg/rcomp/tmp/15p5gk1322582578.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/wessaorg/rcomp/tmp/16q2n21322582578.tab")
+ }
>
> try(system("convert tmp/16qez1322582577.ps tmp/16qez1322582577.png",intern=TRUE))
character(0)
> try(system("convert tmp/28gxb1322582577.ps tmp/28gxb1322582577.png",intern=TRUE))
character(0)
> try(system("convert tmp/33vhf1322582577.ps tmp/33vhf1322582577.png",intern=TRUE))
character(0)
> try(system("convert tmp/4i1eq1322582577.ps tmp/4i1eq1322582577.png",intern=TRUE))
character(0)
> try(system("convert tmp/5zt741322582578.ps tmp/5zt741322582578.png",intern=TRUE))
character(0)
> try(system("convert tmp/6s7ii1322582578.ps tmp/6s7ii1322582578.png",intern=TRUE))
character(0)
> try(system("convert tmp/7pmn91322582578.ps tmp/7pmn91322582578.png",intern=TRUE))
character(0)
> try(system("convert tmp/8vlm21322582578.ps tmp/8vlm21322582578.png",intern=TRUE))
character(0)
> try(system("convert tmp/9b01k1322582578.ps tmp/9b01k1322582578.png",intern=TRUE))
character(0)
> try(system("convert tmp/10hyvk1322582578.ps tmp/10hyvk1322582578.png",intern=TRUE))
character(0)
>
>
> proc.time()
user system elapsed
3.301 0.493 3.949