R version 2.11.1 (2010-05-31)
Copyright (C) 2010 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(1
+ ,28
+ ,6
+ ,6
+ ,6.06
+ ,6.06
+ ,3.53
+ ,3.53
+ ,48
+ ,48
+ ,5
+ ,5
+ ,1
+ ,40
+ ,5
+ ,5
+ ,8.1
+ ,8.1
+ ,4.52
+ ,4.52
+ ,63
+ ,63
+ ,11
+ ,11
+ ,1
+ ,79
+ ,3
+ ,3
+ ,79.38
+ ,79.38
+ ,3.72
+ ,3.72
+ ,113
+ ,113
+ ,13
+ ,13
+ ,1
+ ,16
+ ,2
+ ,2
+ ,26.26
+ ,26.26
+ ,3.17
+ ,3.17
+ ,104
+ ,104
+ ,1
+ ,1
+ ,1
+ ,90
+ ,2
+ ,2
+ ,39.56
+ ,39.56
+ ,3.39
+ ,3.39
+ ,89
+ ,89
+ ,11
+ ,11
+ ,1
+ ,87
+ ,5
+ ,5
+ ,65.61
+ ,65.61
+ ,4.15
+ ,4.15
+ ,97
+ ,97
+ ,3
+ ,3
+ ,1
+ ,53
+ ,5
+ ,5
+ ,80.3
+ ,80.3
+ ,3.09
+ ,3.09
+ ,114
+ ,114
+ ,11
+ ,11
+ ,1
+ ,23
+ ,5
+ ,5
+ ,34.68
+ ,34.68
+ ,2.76
+ ,2.76
+ ,57
+ ,57
+ ,9
+ ,9
+ ,1
+ ,42
+ ,6
+ ,6
+ ,7.17
+ ,7.17
+ ,5.14
+ ,5.14
+ ,127
+ ,127
+ ,10
+ ,10
+ ,1
+ ,64
+ ,4
+ ,4
+ ,65.88
+ ,65.88
+ ,4.78
+ ,4.78
+ ,64
+ ,64
+ ,4
+ ,4
+ ,1
+ ,87
+ ,6
+ ,6
+ ,42.69
+ ,42.69
+ ,4.22
+ ,4.22
+ ,91
+ ,91
+ ,2
+ ,2
+ ,1
+ ,77
+ ,2
+ ,2
+ ,54.94
+ ,54.94
+ ,3.93
+ ,3.93
+ ,127
+ ,127
+ ,2
+ ,2
+ ,1
+ ,70
+ ,4
+ ,4
+ ,89.99
+ ,89.99
+ ,3.01
+ ,3.01
+ ,45
+ ,45
+ ,10
+ ,10
+ ,1
+ ,82
+ ,4
+ ,4
+ ,72.64
+ ,72.64
+ ,5.12
+ ,5.12
+ ,40
+ ,40
+ ,9
+ ,9
+ ,1
+ ,44
+ ,3
+ ,3
+ ,24.96
+ ,24.96
+ ,5.82
+ ,5.82
+ ,33
+ ,33
+ ,1
+ ,1
+ ,1
+ ,36
+ ,2
+ ,2
+ ,57.52
+ ,57.52
+ ,2.83
+ ,2.83
+ ,60
+ ,60
+ ,7
+ ,7
+ ,0
+ ,73
+ ,2
+ ,0
+ ,71.91
+ ,0
+ ,5.11
+ ,0
+ ,50
+ ,0
+ ,3
+ ,0
+ ,0
+ ,75
+ ,3
+ ,0
+ ,65.34
+ ,0
+ ,5.99
+ ,0
+ ,128
+ ,0
+ ,11
+ ,0
+ ,0
+ ,21
+ ,3
+ ,0
+ ,34.62
+ ,0
+ ,3.15
+ ,0
+ ,52
+ ,0
+ ,7
+ ,0
+ ,0
+ ,81
+ ,2
+ ,0
+ ,60.92
+ ,0
+ ,3.5
+ ,0
+ ,40
+ ,0
+ ,1
+ ,0
+ ,0
+ ,99
+ ,3
+ ,0
+ ,56.49
+ ,0
+ ,4.5
+ ,0
+ ,29
+ ,0
+ ,9
+ ,0
+ ,0
+ ,54
+ ,3
+ ,0
+ ,56.19
+ ,0
+ ,3.31
+ ,0
+ ,36
+ ,0
+ ,5
+ ,0
+ ,0
+ ,6
+ ,4
+ ,0
+ ,61.2
+ ,0
+ ,5.31
+ ,0
+ ,49
+ ,0
+ ,9
+ ,0
+ ,0
+ ,71
+ ,5
+ ,0
+ ,58.2
+ ,0
+ ,4.24
+ ,0
+ ,57
+ ,0
+ ,7
+ ,0
+ ,0
+ ,93
+ ,6
+ ,0
+ ,75.91
+ ,0
+ ,5.06
+ ,0
+ ,82
+ ,0
+ ,4
+ ,0
+ ,0
+ ,82
+ ,3
+ ,0
+ ,73.66
+ ,0
+ ,4.72
+ ,0
+ ,34
+ ,0
+ ,10
+ ,0
+ ,0
+ ,32
+ ,4
+ ,0
+ ,73.87
+ ,0
+ ,4.58
+ ,0
+ ,36
+ ,0
+ ,13
+ ,0
+ ,0
+ ,93
+ ,4
+ ,0
+ ,87.21
+ ,0
+ ,5.3
+ ,0
+ ,89
+ ,0
+ ,9
+ ,0
+ ,0
+ ,24
+ ,4
+ ,0
+ ,64.29
+ ,0
+ ,5.11
+ ,0
+ ,69
+ ,0
+ ,5
+ ,0
+ ,0
+ ,96
+ ,5
+ ,0
+ ,71.82
+ ,0
+ ,4.05
+ ,0
+ ,35
+ ,0
+ ,8
+ ,0
+ ,0
+ ,88
+ ,4
+ ,0
+ ,89.31
+ ,0
+ ,4.62
+ ,0
+ ,65
+ ,0
+ ,12
+ ,0
+ ,0
+ ,83
+ ,2
+ ,0
+ ,1.41
+ ,0
+ ,4.66
+ ,0
+ ,70
+ ,0
+ ,8
+ ,0
+ ,0
+ ,23
+ ,6
+ ,0
+ ,35.17
+ ,0
+ ,4.66
+ ,0
+ ,60
+ ,0
+ ,5
+ ,0
+ ,0
+ ,20
+ ,5
+ ,0
+ ,41.08
+ ,0
+ ,5.1
+ ,0
+ ,127
+ ,0
+ ,11
+ ,0
+ ,0
+ ,33
+ ,3
+ ,0
+ ,30.57
+ ,0
+ ,4.97
+ ,0
+ ,96
+ ,0
+ ,8
+ ,0
+ ,0
+ ,88
+ ,2
+ ,0
+ ,68.84
+ ,0
+ ,2.87
+ ,0
+ ,61
+ ,0
+ ,9
+ ,0
+ ,0
+ ,98
+ ,2
+ ,0
+ ,71.05
+ ,0
+ ,4.98
+ ,0
+ ,36
+ ,0
+ ,1
+ ,0
+ ,0
+ ,34
+ ,4
+ ,0
+ ,23.32
+ ,0
+ ,4.55
+ ,0
+ ,55
+ ,0
+ ,9
+ ,0
+ ,0
+ ,59
+ ,3
+ ,0
+ ,61.39
+ ,0
+ ,5.45
+ ,0
+ ,75
+ ,0
+ ,2
+ ,0
+ ,0
+ ,26
+ ,6
+ ,0
+ ,8.41
+ ,0
+ ,4.36
+ ,0
+ ,42
+ ,0
+ ,3
+ ,0
+ ,0
+ ,13
+ ,1
+ ,0
+ ,64.06
+ ,0
+ ,4.74
+ ,0
+ ,83
+ ,0
+ ,3
+ ,0
+ ,0
+ ,6
+ ,2
+ ,0
+ ,26.8
+ ,0
+ ,5.44
+ ,0
+ ,56
+ ,0
+ ,1
+ ,0
+ ,0
+ ,49
+ ,4
+ ,0
+ ,12.78
+ ,0
+ ,5.78
+ ,0
+ ,114
+ ,0
+ ,5
+ ,0
+ ,0
+ ,3
+ ,5
+ ,0
+ ,23.84
+ ,0
+ ,2.92
+ ,0
+ ,33
+ ,0
+ ,4
+ ,0
+ ,0
+ ,76
+ ,4
+ ,0
+ ,5.68
+ ,0
+ ,3.22
+ ,0
+ ,80
+ ,0
+ ,6
+ ,0
+ ,0
+ ,12
+ ,2
+ ,0
+ ,45.92
+ ,0
+ ,3.04
+ ,0
+ ,115
+ ,0
+ ,7
+ ,0
+ ,0
+ ,63
+ ,5
+ ,0
+ ,18.17
+ ,0
+ ,3.11
+ ,0
+ ,127
+ ,0
+ ,13
+ ,0
+ ,0
+ ,35
+ ,1
+ ,0
+ ,29.12
+ ,0
+ ,3.87
+ ,0
+ ,45
+ ,0
+ ,9
+ ,0
+ ,0
+ ,69
+ ,1
+ ,0
+ ,40.08
+ ,0
+ ,3.75
+ ,0
+ ,74
+ ,0
+ ,11
+ ,0
+ ,0
+ ,10
+ ,5
+ ,0
+ ,1.08
+ ,0
+ ,4.82
+ ,0
+ ,105
+ ,0
+ ,10
+ ,0)
+ ,dim=c(12
+ ,50)
+ ,dimnames=list(c('klant'
+ ,'slaagkans'
+ ,'verzekeraar'
+ ,'verzekeraar_klant'
+ ,'kost'
+ ,'kost_klant'
+ ,'grootte'
+ ,'grootte_klant'
+ ,'snelheid'
+ ,'snelheid_klant'
+ ,'maand'
+ ,'maand_klant')
+ ,1:50))
> y <- array(NA,dim=c(12,50),dimnames=list(c('klant','slaagkans','verzekeraar','verzekeraar_klant','kost','kost_klant','grootte','grootte_klant','snelheid','snelheid_klant','maand','maand_klant'),1:50))
> 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 = '2'
> #'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
slaagkans klant verzekeraar verzekeraar_klant kost kost_klant grootte
1 28 1 6 6 6.06 6.06 3.53
2 40 1 5 5 8.10 8.10 4.52
3 79 1 3 3 79.38 79.38 3.72
4 16 1 2 2 26.26 26.26 3.17
5 90 1 2 2 39.56 39.56 3.39
6 87 1 5 5 65.61 65.61 4.15
7 53 1 5 5 80.30 80.30 3.09
8 23 1 5 5 34.68 34.68 2.76
9 42 1 6 6 7.17 7.17 5.14
10 64 1 4 4 65.88 65.88 4.78
11 87 1 6 6 42.69 42.69 4.22
12 77 1 2 2 54.94 54.94 3.93
13 70 1 4 4 89.99 89.99 3.01
14 82 1 4 4 72.64 72.64 5.12
15 44 1 3 3 24.96 24.96 5.82
16 36 1 2 2 57.52 57.52 2.83
17 73 0 2 0 71.91 0.00 5.11
18 75 0 3 0 65.34 0.00 5.99
19 21 0 3 0 34.62 0.00 3.15
20 81 0 2 0 60.92 0.00 3.50
21 99 0 3 0 56.49 0.00 4.50
22 54 0 3 0 56.19 0.00 3.31
23 6 0 4 0 61.20 0.00 5.31
24 71 0 5 0 58.20 0.00 4.24
25 93 0 6 0 75.91 0.00 5.06
26 82 0 3 0 73.66 0.00 4.72
27 32 0 4 0 73.87 0.00 4.58
28 93 0 4 0 87.21 0.00 5.30
29 24 0 4 0 64.29 0.00 5.11
30 96 0 5 0 71.82 0.00 4.05
31 88 0 4 0 89.31 0.00 4.62
32 83 0 2 0 1.41 0.00 4.66
33 23 0 6 0 35.17 0.00 4.66
34 20 0 5 0 41.08 0.00 5.10
35 33 0 3 0 30.57 0.00 4.97
36 88 0 2 0 68.84 0.00 2.87
37 98 0 2 0 71.05 0.00 4.98
38 34 0 4 0 23.32 0.00 4.55
39 59 0 3 0 61.39 0.00 5.45
40 26 0 6 0 8.41 0.00 4.36
41 13 0 1 0 64.06 0.00 4.74
42 6 0 2 0 26.80 0.00 5.44
43 49 0 4 0 12.78 0.00 5.78
44 3 0 5 0 23.84 0.00 2.92
45 76 0 4 0 5.68 0.00 3.22
46 12 0 2 0 45.92 0.00 3.04
47 63 0 5 0 18.17 0.00 3.11
48 35 0 1 0 29.12 0.00 3.87
49 69 0 1 0 40.08 0.00 3.75
50 10 0 5 0 1.08 0.00 4.82
grootte_klant snelheid snelheid_klant maand maand_klant
1 3.53 48 48 5 5
2 4.52 63 63 11 11
3 3.72 113 113 13 13
4 3.17 104 104 1 1
5 3.39 89 89 11 11
6 4.15 97 97 3 3
7 3.09 114 114 11 11
8 2.76 57 57 9 9
9 5.14 127 127 10 10
10 4.78 64 64 4 4
11 4.22 91 91 2 2
12 3.93 127 127 2 2
13 3.01 45 45 10 10
14 5.12 40 40 9 9
15 5.82 33 33 1 1
16 2.83 60 60 7 7
17 0.00 50 0 3 0
18 0.00 128 0 11 0
19 0.00 52 0 7 0
20 0.00 40 0 1 0
21 0.00 29 0 9 0
22 0.00 36 0 5 0
23 0.00 49 0 9 0
24 0.00 57 0 7 0
25 0.00 82 0 4 0
26 0.00 34 0 10 0
27 0.00 36 0 13 0
28 0.00 89 0 9 0
29 0.00 69 0 5 0
30 0.00 35 0 8 0
31 0.00 65 0 12 0
32 0.00 70 0 8 0
33 0.00 60 0 5 0
34 0.00 127 0 11 0
35 0.00 96 0 8 0
36 0.00 61 0 9 0
37 0.00 36 0 1 0
38 0.00 55 0 9 0
39 0.00 75 0 2 0
40 0.00 42 0 3 0
41 0.00 83 0 3 0
42 0.00 56 0 1 0
43 0.00 114 0 5 0
44 0.00 33 0 4 0
45 0.00 80 0 6 0
46 0.00 115 0 7 0
47 0.00 127 0 13 0
48 0.00 45 0 9 0
49 0.00 74 0 11 0
50 0.00 105 0 10 0
> k <- length(x[1,])
> df <- as.data.frame(x)
> (mylm <- lm(df))
Call:
lm(formula = df)
Coefficients:
(Intercept) klant verzekeraar verzekeraar_klant
36.78638 -70.40610 -1.04099 1.39545
kost kost_klant grootte grootte_klant
0.55479 0.05675 -2.24036 13.59903
snelheid snelheid_klant maand maand_klant
-0.04266 0.22118 0.85427 -0.61258
> (mysum <- summary(mylm))
Call:
lm(formula = df)
Residuals:
Min 1Q Median 3Q Max
-54.277 -15.175 1.977 15.978 54.106
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 36.78638 30.30401 1.214 0.2323
klant -70.40610 59.29137 -1.187 0.2424
verzekeraar -1.04099 3.52878 -0.295 0.7696
verzekeraar_klant 1.39545 6.37253 0.219 0.8278
kost 0.55479 0.21581 2.571 0.0142 *
kost_klant 0.05675 0.36290 0.156 0.8766
grootte -2.24036 6.31464 -0.355 0.7247
grootte_klant 13.59903 10.75192 1.265 0.2136
snelheid -0.04266 0.19360 -0.220 0.8268
snelheid_klant 0.22118 0.30496 0.725 0.4727
maand 0.85427 1.55729 0.549 0.5865
maand_klant -0.61258 2.47010 -0.248 0.8055
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 28.8 on 38 degrees of freedom
Multiple R-squared: 0.2858, Adjusted R-squared: 0.07908
F-statistic: 1.383 on 11 and 38 DF, p-value: 0.2208
> 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.67388239 0.65223522 0.3261176
[2,] 0.52253176 0.95493649 0.4774682
[3,] 0.36672276 0.73344552 0.6332772
[4,] 0.24138290 0.48276579 0.7586171
[5,] 0.15692286 0.31384573 0.8430771
[6,] 0.09239675 0.18479350 0.9076033
[7,] 0.05889686 0.11779372 0.9411031
[8,] 0.02992166 0.05984332 0.9700783
[9,] 0.02847258 0.05694515 0.9715274
[10,] 0.09671493 0.19342987 0.9032851
[11,] 0.08325173 0.16650346 0.9167483
[12,] 0.05862842 0.11725685 0.9413716
[13,] 0.09677436 0.19354872 0.9032256
[14,] 0.06433422 0.12866844 0.9356658
[15,] 0.06728468 0.13456935 0.9327153
[16,] 0.05195056 0.10390113 0.9480494
[17,] 0.03023170 0.06046340 0.9697683
[18,] 0.06943095 0.13886189 0.9305691
[19,] 0.04222474 0.08444948 0.9577753
[20,] 0.03432928 0.06865855 0.9656707
[21,] 0.01537506 0.03075012 0.9846249
> postscript(file="/var/www/rcomp/tmp/117yk1290537364.ps",horizontal=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/2ugx51290537364.ps",horizontal=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/3ugx51290537364.ps",horizontal=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/4ugx51290537364.ps",horizontal=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/5ugx51290537364.ps",horizontal=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 = 50
Frequency = 1
1 2 3 4 5 6
5.9135175 1.6474533 -2.5569342 -21.9630338 41.6654450 13.5440374
7 8 9 10 11 12
-22.3676766 -10.0616728 -14.3641874 -10.7731294 27.7238761 8.5176422
13 14 15 16 17 18
2.5291144 2.3069113 -10.9481076 -10.8132555 9.4189773 14.5699550
19 20 21 22 23 24
-28.5746612 21.1910856 37.6267543 -6.1571387 -54.2774071 13.0805990
25 26 27 28 29 30
31.7626428 10.9528846 -40.9137611 19.9765000 -34.1694889 28.3058289
31 32 33 34 35 36
8.7013103 54.1056078 -18.3240694 -26.9253825 -9.2274455 16.4474678
37 38 39 40 41 42
35.7161202 -6.7085608 4.9789152 -0.2093026 -46.6879812 -29.8505007
43 44 45 46 47 48
22.8287048 -37.2750804 45.7276473 -42.4435547 22.6180714 -13.9993792
49 50
13.1799076 -15.4452658
> postscript(file="/var/www/rcomp/tmp/6npw81290537364.ps",horizontal=F,pagecentre=F,paper="special",width=8.3333333333333,height=5.5555555555556)
> dum <- cbind(lag(myerror,k=1),myerror)
> dum
Time Series:
Start = 0
End = 50
Frequency = 1
lag(myerror, k = 1) myerror
0 5.9135175 NA
1 1.6474533 5.9135175
2 -2.5569342 1.6474533
3 -21.9630338 -2.5569342
4 41.6654450 -21.9630338
5 13.5440374 41.6654450
6 -22.3676766 13.5440374
7 -10.0616728 -22.3676766
8 -14.3641874 -10.0616728
9 -10.7731294 -14.3641874
10 27.7238761 -10.7731294
11 8.5176422 27.7238761
12 2.5291144 8.5176422
13 2.3069113 2.5291144
14 -10.9481076 2.3069113
15 -10.8132555 -10.9481076
16 9.4189773 -10.8132555
17 14.5699550 9.4189773
18 -28.5746612 14.5699550
19 21.1910856 -28.5746612
20 37.6267543 21.1910856
21 -6.1571387 37.6267543
22 -54.2774071 -6.1571387
23 13.0805990 -54.2774071
24 31.7626428 13.0805990
25 10.9528846 31.7626428
26 -40.9137611 10.9528846
27 19.9765000 -40.9137611
28 -34.1694889 19.9765000
29 28.3058289 -34.1694889
30 8.7013103 28.3058289
31 54.1056078 8.7013103
32 -18.3240694 54.1056078
33 -26.9253825 -18.3240694
34 -9.2274455 -26.9253825
35 16.4474678 -9.2274455
36 35.7161202 16.4474678
37 -6.7085608 35.7161202
38 4.9789152 -6.7085608
39 -0.2093026 4.9789152
40 -46.6879812 -0.2093026
41 -29.8505007 -46.6879812
42 22.8287048 -29.8505007
43 -37.2750804 22.8287048
44 45.7276473 -37.2750804
45 -42.4435547 45.7276473
46 22.6180714 -42.4435547
47 -13.9993792 22.6180714
48 13.1799076 -13.9993792
49 -15.4452658 13.1799076
50 NA -15.4452658
> dum1 <- dum[2:length(myerror),]
> dum1
lag(myerror, k = 1) myerror
[1,] 1.6474533 5.9135175
[2,] -2.5569342 1.6474533
[3,] -21.9630338 -2.5569342
[4,] 41.6654450 -21.9630338
[5,] 13.5440374 41.6654450
[6,] -22.3676766 13.5440374
[7,] -10.0616728 -22.3676766
[8,] -14.3641874 -10.0616728
[9,] -10.7731294 -14.3641874
[10,] 27.7238761 -10.7731294
[11,] 8.5176422 27.7238761
[12,] 2.5291144 8.5176422
[13,] 2.3069113 2.5291144
[14,] -10.9481076 2.3069113
[15,] -10.8132555 -10.9481076
[16,] 9.4189773 -10.8132555
[17,] 14.5699550 9.4189773
[18,] -28.5746612 14.5699550
[19,] 21.1910856 -28.5746612
[20,] 37.6267543 21.1910856
[21,] -6.1571387 37.6267543
[22,] -54.2774071 -6.1571387
[23,] 13.0805990 -54.2774071
[24,] 31.7626428 13.0805990
[25,] 10.9528846 31.7626428
[26,] -40.9137611 10.9528846
[27,] 19.9765000 -40.9137611
[28,] -34.1694889 19.9765000
[29,] 28.3058289 -34.1694889
[30,] 8.7013103 28.3058289
[31,] 54.1056078 8.7013103
[32,] -18.3240694 54.1056078
[33,] -26.9253825 -18.3240694
[34,] -9.2274455 -26.9253825
[35,] 16.4474678 -9.2274455
[36,] 35.7161202 16.4474678
[37,] -6.7085608 35.7161202
[38,] 4.9789152 -6.7085608
[39,] -0.2093026 4.9789152
[40,] -46.6879812 -0.2093026
[41,] -29.8505007 -46.6879812
[42,] 22.8287048 -29.8505007
[43,] -37.2750804 22.8287048
[44,] 45.7276473 -37.2750804
[45,] -42.4435547 45.7276473
[46,] 22.6180714 -42.4435547
[47,] -13.9993792 22.6180714
[48,] 13.1799076 -13.9993792
[49,] -15.4452658 13.1799076
> z <- as.data.frame(dum1)
> z
lag(myerror, k = 1) myerror
1 1.6474533 5.9135175
2 -2.5569342 1.6474533
3 -21.9630338 -2.5569342
4 41.6654450 -21.9630338
5 13.5440374 41.6654450
6 -22.3676766 13.5440374
7 -10.0616728 -22.3676766
8 -14.3641874 -10.0616728
9 -10.7731294 -14.3641874
10 27.7238761 -10.7731294
11 8.5176422 27.7238761
12 2.5291144 8.5176422
13 2.3069113 2.5291144
14 -10.9481076 2.3069113
15 -10.8132555 -10.9481076
16 9.4189773 -10.8132555
17 14.5699550 9.4189773
18 -28.5746612 14.5699550
19 21.1910856 -28.5746612
20 37.6267543 21.1910856
21 -6.1571387 37.6267543
22 -54.2774071 -6.1571387
23 13.0805990 -54.2774071
24 31.7626428 13.0805990
25 10.9528846 31.7626428
26 -40.9137611 10.9528846
27 19.9765000 -40.9137611
28 -34.1694889 19.9765000
29 28.3058289 -34.1694889
30 8.7013103 28.3058289
31 54.1056078 8.7013103
32 -18.3240694 54.1056078
33 -26.9253825 -18.3240694
34 -9.2274455 -26.9253825
35 16.4474678 -9.2274455
36 35.7161202 16.4474678
37 -6.7085608 35.7161202
38 4.9789152 -6.7085608
39 -0.2093026 4.9789152
40 -46.6879812 -0.2093026
41 -29.8505007 -46.6879812
42 22.8287048 -29.8505007
43 -37.2750804 22.8287048
44 45.7276473 -37.2750804
45 -42.4435547 45.7276473
46 22.6180714 -42.4435547
47 -13.9993792 22.6180714
48 13.1799076 -13.9993792
49 -15.4452658 13.1799076
> 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/7gzdb1290537364.ps",horizontal=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/8gzdb1290537364.ps",horizontal=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/9gzdb1290537364.ps",horizontal=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/1088ve1290537364.ps",horizontal=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/11uqt21290537364.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/12f9sq1290537364.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/13mspj1290537364.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/14fj651290537364.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/150kms1290537364.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/16wbkj1290537364.tab")
+ }
>
> try(system("convert tmp/117yk1290537364.ps tmp/117yk1290537364.png",intern=TRUE))
character(0)
> try(system("convert tmp/2ugx51290537364.ps tmp/2ugx51290537364.png",intern=TRUE))
character(0)
> try(system("convert tmp/3ugx51290537364.ps tmp/3ugx51290537364.png",intern=TRUE))
character(0)
> try(system("convert tmp/4ugx51290537364.ps tmp/4ugx51290537364.png",intern=TRUE))
character(0)
> try(system("convert tmp/5ugx51290537364.ps tmp/5ugx51290537364.png",intern=TRUE))
character(0)
> try(system("convert tmp/6npw81290537364.ps tmp/6npw81290537364.png",intern=TRUE))
character(0)
> try(system("convert tmp/7gzdb1290537364.ps tmp/7gzdb1290537364.png",intern=TRUE))
character(0)
> try(system("convert tmp/8gzdb1290537364.ps tmp/8gzdb1290537364.png",intern=TRUE))
character(0)
> try(system("convert tmp/9gzdb1290537364.ps tmp/9gzdb1290537364.png",intern=TRUE))
character(0)
> try(system("convert tmp/1088ve1290537364.ps tmp/1088ve1290537364.png",intern=TRUE))
character(0)
>
>
> proc.time()
user system elapsed
3.580 2.090 5.634