R version 2.8.0 (2008-10-20)
Copyright (C) 2008 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.
Natural language support but running in an English locale
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(73
+ ,2
+ ,71.91
+ ,5.11
+ ,50
+ ,3
+ ,28
+ ,6
+ ,6.06
+ ,3.53
+ ,48
+ ,5
+ ,40
+ ,5
+ ,8.1
+ ,4.52
+ ,63
+ ,11
+ ,79
+ ,3
+ ,79.38
+ ,3.72
+ ,113
+ ,13
+ ,75
+ ,3
+ ,65.34
+ ,5.99
+ ,128
+ ,11
+ ,21
+ ,3
+ ,34.62
+ ,3.15
+ ,52
+ ,7
+ ,16
+ ,2
+ ,26.26
+ ,3.17
+ ,104
+ ,1
+ ,81
+ ,2
+ ,60.92
+ ,3.5
+ ,40
+ ,1
+ ,90
+ ,2
+ ,39.56
+ ,3.39
+ ,89
+ ,11
+ ,87
+ ,5
+ ,65.61
+ ,4.15
+ ,97
+ ,3
+ ,99
+ ,3
+ ,56.49
+ ,4.5
+ ,29
+ ,9
+ ,54
+ ,3
+ ,56.19
+ ,3.31
+ ,36
+ ,5
+ ,53
+ ,5
+ ,80.3
+ ,3.09
+ ,114
+ ,11
+ ,6
+ ,4
+ ,61.2
+ ,5.31
+ ,49
+ ,9
+ ,71
+ ,5
+ ,58.2
+ ,4.24
+ ,57
+ ,7
+ ,93
+ ,6
+ ,75.91
+ ,5.06
+ ,82
+ ,4
+ ,82
+ ,3
+ ,73.66
+ ,4.72
+ ,34
+ ,10
+ ,32
+ ,4
+ ,73.87
+ ,4.58
+ ,36
+ ,13
+ ,93
+ ,4
+ ,87.21
+ ,5.3
+ ,89
+ ,9
+ ,24
+ ,4
+ ,64.29
+ ,5.11
+ ,69
+ ,5
+ ,96
+ ,5
+ ,71.82
+ ,4.05
+ ,35
+ ,8
+ ,88
+ ,4
+ ,89.31
+ ,4.62
+ ,65
+ ,12
+ ,83
+ ,2
+ ,1.41
+ ,4.66
+ ,70
+ ,8
+ ,23
+ ,6
+ ,35.17
+ ,4.66
+ ,60
+ ,5
+ ,23
+ ,5
+ ,34.68
+ ,2.76
+ ,57
+ ,9
+ ,20
+ ,5
+ ,41.08
+ ,5.1
+ ,127
+ ,11
+ ,33
+ ,3
+ ,30.57
+ ,4.97
+ ,96
+ ,8
+ ,88
+ ,2
+ ,68.84
+ ,2.87
+ ,61
+ ,9
+ ,42
+ ,6
+ ,7.17
+ ,5.14
+ ,127
+ ,10
+ ,98
+ ,2
+ ,71.05
+ ,4.98
+ ,36
+ ,1
+ ,34
+ ,4
+ ,23.32
+ ,4.55
+ ,55
+ ,9
+ ,59
+ ,3
+ ,61.39
+ ,5.45
+ ,75
+ ,2
+ ,26
+ ,6
+ ,8.41
+ ,4.36
+ ,42
+ ,3
+ ,64
+ ,4
+ ,65.88
+ ,4.78
+ ,64
+ ,4
+ ,13
+ ,1
+ ,64.06
+ ,4.74
+ ,83
+ ,3
+ ,6
+ ,2
+ ,26.8
+ ,5.44
+ ,56
+ ,1
+ ,49
+ ,4
+ ,12.78
+ ,5.78
+ ,114
+ ,5
+ ,3
+ ,5
+ ,23.84
+ ,2.92
+ ,33
+ ,4
+ ,87
+ ,6
+ ,42.69
+ ,4.22
+ ,91
+ ,2
+ ,77
+ ,2
+ ,54.94
+ ,3.93
+ ,127
+ ,2
+ ,70
+ ,4
+ ,89.99
+ ,3.01
+ ,45
+ ,10
+ ,76
+ ,4
+ ,5.68
+ ,3.22
+ ,80
+ ,6
+ ,82
+ ,4
+ ,72.64
+ ,5.12
+ ,40
+ ,9
+ ,12
+ ,2
+ ,45.92
+ ,3.04
+ ,115
+ ,7
+ ,44
+ ,3
+ ,24.96
+ ,5.82
+ ,33
+ ,1
+ ,63
+ ,5
+ ,18.17
+ ,3.11
+ ,127
+ ,13
+ ,35
+ ,1
+ ,29.12
+ ,3.87
+ ,45
+ ,9
+ ,69
+ ,1
+ ,40.08
+ ,3.75
+ ,74
+ ,11
+ ,10
+ ,5
+ ,1.08
+ ,4.82
+ ,105
+ ,10
+ ,36
+ ,2
+ ,57.52
+ ,2.83
+ ,60
+ ,7)
+ ,dim=c(6
+ ,50)
+ ,dimnames=list(c('slaagkans'
+ ,'verzekeraar'
+ ,'kost'
+ ,'grootte'
+ ,'snelheid'
+ ,'maand')
+ ,1:50))
> y <- array(NA,dim=c(6,50),dimnames=list(c('slaagkans','verzekeraar','kost','grootte','snelheid','maand'),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 = '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
slaagkans verzekeraar kost grootte snelheid maand t
1 73 2 71.91 5.11 50 3 1
2 28 6 6.06 3.53 48 5 2
3 40 5 8.10 4.52 63 11 3
4 79 3 79.38 3.72 113 13 4
5 75 3 65.34 5.99 128 11 5
6 21 3 34.62 3.15 52 7 6
7 16 2 26.26 3.17 104 1 7
8 81 2 60.92 3.50 40 1 8
9 90 2 39.56 3.39 89 11 9
10 87 5 65.61 4.15 97 3 10
11 99 3 56.49 4.50 29 9 11
12 54 3 56.19 3.31 36 5 12
13 53 5 80.30 3.09 114 11 13
14 6 4 61.20 5.31 49 9 14
15 71 5 58.20 4.24 57 7 15
16 93 6 75.91 5.06 82 4 16
17 82 3 73.66 4.72 34 10 17
18 32 4 73.87 4.58 36 13 18
19 93 4 87.21 5.30 89 9 19
20 24 4 64.29 5.11 69 5 20
21 96 5 71.82 4.05 35 8 21
22 88 4 89.31 4.62 65 12 22
23 83 2 1.41 4.66 70 8 23
24 23 6 35.17 4.66 60 5 24
25 23 5 34.68 2.76 57 9 25
26 20 5 41.08 5.10 127 11 26
27 33 3 30.57 4.97 96 8 27
28 88 2 68.84 2.87 61 9 28
29 42 6 7.17 5.14 127 10 29
30 98 2 71.05 4.98 36 1 30
31 34 4 23.32 4.55 55 9 31
32 59 3 61.39 5.45 75 2 32
33 26 6 8.41 4.36 42 3 33
34 64 4 65.88 4.78 64 4 34
35 13 1 64.06 4.74 83 3 35
36 6 2 26.80 5.44 56 1 36
37 49 4 12.78 5.78 114 5 37
38 3 5 23.84 2.92 33 4 38
39 87 6 42.69 4.22 91 2 39
40 77 2 54.94 3.93 127 2 40
41 70 4 89.99 3.01 45 10 41
42 76 4 5.68 3.22 80 6 42
43 82 4 72.64 5.12 40 9 43
44 12 2 45.92 3.04 115 7 44
45 44 3 24.96 5.82 33 1 45
46 63 5 18.17 3.11 127 13 46
47 35 1 29.12 3.87 45 9 47
48 69 1 40.08 3.75 74 11 48
49 10 5 1.08 4.82 105 10 49
50 36 2 57.52 2.83 60 7 50
> k <- length(x[1,])
> df <- as.data.frame(x)
> (mylm <- lm(df))
Call:
lm(formula = df)
Coefficients:
(Intercept) verzekeraar kost grootte snelheid maand
26.83008 -0.39719 0.52564 0.73477 0.03087 0.31850
t
-0.14492
> (mysum <- summary(mylm))
Call:
lm(formula = df)
Residuals:
Min 1Q Median 3Q Max
-57.662 -20.893 2.429 19.961 51.423
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 26.83008 27.27626 0.984 0.33079
verzekeraar -0.39719 2.82093 -0.141 0.88869
kost 0.52564 0.16561 3.174 0.00278 **
grootte 0.73477 4.48739 0.164 0.87070
snelheid 0.03087 0.13552 0.228 0.82090
maand 0.31850 1.14555 0.278 0.78232
t -0.14492 0.28753 -0.504 0.61682
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Residual standard error: 28.03 on 43 degrees of freedom
Multiple R-squared: 0.2347, Adjusted R-squared: 0.1279
F-statistic: 2.198 on 6 and 43 DF, p-value: 0.06156
> 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.4093732 0.8187465 0.5906268
[2,] 0.2966289 0.5932579 0.7033711
[3,] 0.3187313 0.6374626 0.6812687
[4,] 0.3192010 0.6384019 0.6807990
[5,] 0.7557064 0.4885871 0.2442936
[6,] 0.6923334 0.6153331 0.3076666
[7,] 0.6589178 0.6821645 0.3410822
[8,] 0.5653993 0.8692014 0.4346007
[9,] 0.6088831 0.7822337 0.3911169
[10,] 0.5324224 0.9351552 0.4675776
[11,] 0.5680940 0.8638119 0.4319060
[12,] 0.5741255 0.8517491 0.4258745
[13,] 0.4910720 0.9821441 0.5089280
[14,] 0.6858278 0.6283443 0.3141722
[15,] 0.6535389 0.6929222 0.3464611
[16,] 0.6032726 0.7934549 0.3967274
[17,] 0.6313581 0.7372838 0.3686419
[18,] 0.5570707 0.8858585 0.4429293
[19,] 0.5272896 0.9454208 0.4727104
[20,] 0.4527330 0.9054661 0.5472670
[21,] 0.5495178 0.9009644 0.4504822
[22,] 0.4543791 0.9087581 0.5456209
[23,] 0.3647723 0.7295447 0.6352277
[24,] 0.2719279 0.5438558 0.7280721
[25,] 0.1934681 0.3869361 0.8065319
[26,] 0.2593985 0.5187969 0.7406015
[27,] 0.2764425 0.5528851 0.7235575
[28,] 0.2338183 0.4676367 0.7661817
[29,] 0.4471361 0.8942723 0.5528639
[30,] 0.3749197 0.7498394 0.6250803
[31,] 0.3760145 0.7520291 0.6239855
> postscript(file="/var/www/html/freestat/rcomp/tmp/1w4hg1290537549.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/html/freestat/rcomp/tmp/2w4hg1290537549.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/html/freestat/rcomp/tmp/3w4hg1290537549.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/html/freestat/rcomp/tmp/4ovyj1290537549.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/html/freestat/rcomp/tmp/5ovyj1290537549.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 7
3.057196 -5.010396 2.563579 1.854224 3.885120 -28.115691 -28.682494
8 9 10 11 12 13 14
19.977023 35.732750 22.119081 38.194299 -4.570791 -21.461678 -57.662031
15 16 17 18 19 20 21
10.633236 23.447627 12.404161 -38.078494 15.163396 -39.613139 29.843809
22 23 24 25 26 27 28
9.779297 51.423456 -23.324130 -23.104186 -33.840489 -13.957565 22.979012
29 30 31 32 33 34 35
7.104857 33.876574 -6.914137 -1.226481 -3.540793 2.295283 -49.033312
36 37 38 39 40 41 42
-34.949917 13.044574 -33.306560 39.218755 20.437709 -3.387327 47.113208
43 44 45 46 47 48 49
16.944772 -39.809534 6.149474 24.925370 -7.027354 19.912581 -18.278601
50
-19.185318
> postscript(file="/var/www/html/freestat/rcomp/tmp/6hmfm1290537549.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 3.057196 NA
1 -5.010396 3.057196
2 2.563579 -5.010396
3 1.854224 2.563579
4 3.885120 1.854224
5 -28.115691 3.885120
6 -28.682494 -28.115691
7 19.977023 -28.682494
8 35.732750 19.977023
9 22.119081 35.732750
10 38.194299 22.119081
11 -4.570791 38.194299
12 -21.461678 -4.570791
13 -57.662031 -21.461678
14 10.633236 -57.662031
15 23.447627 10.633236
16 12.404161 23.447627
17 -38.078494 12.404161
18 15.163396 -38.078494
19 -39.613139 15.163396
20 29.843809 -39.613139
21 9.779297 29.843809
22 51.423456 9.779297
23 -23.324130 51.423456
24 -23.104186 -23.324130
25 -33.840489 -23.104186
26 -13.957565 -33.840489
27 22.979012 -13.957565
28 7.104857 22.979012
29 33.876574 7.104857
30 -6.914137 33.876574
31 -1.226481 -6.914137
32 -3.540793 -1.226481
33 2.295283 -3.540793
34 -49.033312 2.295283
35 -34.949917 -49.033312
36 13.044574 -34.949917
37 -33.306560 13.044574
38 39.218755 -33.306560
39 20.437709 39.218755
40 -3.387327 20.437709
41 47.113208 -3.387327
42 16.944772 47.113208
43 -39.809534 16.944772
44 6.149474 -39.809534
45 24.925370 6.149474
46 -7.027354 24.925370
47 19.912581 -7.027354
48 -18.278601 19.912581
49 -19.185318 -18.278601
50 NA -19.185318
> dum1 <- dum[2:length(myerror),]
> dum1
lag(myerror, k = 1) myerror
[1,] -5.010396 3.057196
[2,] 2.563579 -5.010396
[3,] 1.854224 2.563579
[4,] 3.885120 1.854224
[5,] -28.115691 3.885120
[6,] -28.682494 -28.115691
[7,] 19.977023 -28.682494
[8,] 35.732750 19.977023
[9,] 22.119081 35.732750
[10,] 38.194299 22.119081
[11,] -4.570791 38.194299
[12,] -21.461678 -4.570791
[13,] -57.662031 -21.461678
[14,] 10.633236 -57.662031
[15,] 23.447627 10.633236
[16,] 12.404161 23.447627
[17,] -38.078494 12.404161
[18,] 15.163396 -38.078494
[19,] -39.613139 15.163396
[20,] 29.843809 -39.613139
[21,] 9.779297 29.843809
[22,] 51.423456 9.779297
[23,] -23.324130 51.423456
[24,] -23.104186 -23.324130
[25,] -33.840489 -23.104186
[26,] -13.957565 -33.840489
[27,] 22.979012 -13.957565
[28,] 7.104857 22.979012
[29,] 33.876574 7.104857
[30,] -6.914137 33.876574
[31,] -1.226481 -6.914137
[32,] -3.540793 -1.226481
[33,] 2.295283 -3.540793
[34,] -49.033312 2.295283
[35,] -34.949917 -49.033312
[36,] 13.044574 -34.949917
[37,] -33.306560 13.044574
[38,] 39.218755 -33.306560
[39,] 20.437709 39.218755
[40,] -3.387327 20.437709
[41,] 47.113208 -3.387327
[42,] 16.944772 47.113208
[43,] -39.809534 16.944772
[44,] 6.149474 -39.809534
[45,] 24.925370 6.149474
[46,] -7.027354 24.925370
[47,] 19.912581 -7.027354
[48,] -18.278601 19.912581
[49,] -19.185318 -18.278601
> z <- as.data.frame(dum1)
> z
lag(myerror, k = 1) myerror
1 -5.010396 3.057196
2 2.563579 -5.010396
3 1.854224 2.563579
4 3.885120 1.854224
5 -28.115691 3.885120
6 -28.682494 -28.115691
7 19.977023 -28.682494
8 35.732750 19.977023
9 22.119081 35.732750
10 38.194299 22.119081
11 -4.570791 38.194299
12 -21.461678 -4.570791
13 -57.662031 -21.461678
14 10.633236 -57.662031
15 23.447627 10.633236
16 12.404161 23.447627
17 -38.078494 12.404161
18 15.163396 -38.078494
19 -39.613139 15.163396
20 29.843809 -39.613139
21 9.779297 29.843809
22 51.423456 9.779297
23 -23.324130 51.423456
24 -23.104186 -23.324130
25 -33.840489 -23.104186
26 -13.957565 -33.840489
27 22.979012 -13.957565
28 7.104857 22.979012
29 33.876574 7.104857
30 -6.914137 33.876574
31 -1.226481 -6.914137
32 -3.540793 -1.226481
33 2.295283 -3.540793
34 -49.033312 2.295283
35 -34.949917 -49.033312
36 13.044574 -34.949917
37 -33.306560 13.044574
38 39.218755 -33.306560
39 20.437709 39.218755
40 -3.387327 20.437709
41 47.113208 -3.387327
42 16.944772 47.113208
43 -39.809534 16.944772
44 6.149474 -39.809534
45 24.925370 6.149474
46 -7.027354 24.925370
47 19.912581 -7.027354
48 -18.278601 19.912581
49 -19.185318 -18.278601
> 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/freestat/rcomp/tmp/7hmfm1290537549.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/html/freestat/rcomp/tmp/8sde71290537549.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/html/freestat/rcomp/tmp/9sde71290537549.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/html/freestat/rcomp/tmp/102nea1290537549.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/html/freestat/rcomp/createtable file can be downloaded at http://www.wessa.net/cretab
> load(file="/var/www/html/freestat/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/freestat/rcomp/tmp/116nuf1290537549.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/freestat/rcomp/tmp/129ob31290537549.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/freestat/rcomp/tmp/13ny8c1290537549.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/freestat/rcomp/tmp/149y7i1290537549.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/freestat/rcomp/tmp/15chno1290537549.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/freestat/rcomp/tmp/168qlf1290537549.tab")
+ }
>
> try(system("convert tmp/1w4hg1290537549.ps tmp/1w4hg1290537549.png",intern=TRUE))
character(0)
> try(system("convert tmp/2w4hg1290537549.ps tmp/2w4hg1290537549.png",intern=TRUE))
character(0)
> try(system("convert tmp/3w4hg1290537549.ps tmp/3w4hg1290537549.png",intern=TRUE))
character(0)
> try(system("convert tmp/4ovyj1290537549.ps tmp/4ovyj1290537549.png",intern=TRUE))
character(0)
> try(system("convert tmp/5ovyj1290537549.ps tmp/5ovyj1290537549.png",intern=TRUE))
character(0)
> try(system("convert tmp/6hmfm1290537549.ps tmp/6hmfm1290537549.png",intern=TRUE))
character(0)
> try(system("convert tmp/7hmfm1290537549.ps tmp/7hmfm1290537549.png",intern=TRUE))
character(0)
> try(system("convert tmp/8sde71290537549.ps tmp/8sde71290537549.png",intern=TRUE))
character(0)
> try(system("convert tmp/9sde71290537549.ps tmp/9sde71290537549.png",intern=TRUE))
character(0)
> try(system("convert tmp/102nea1290537549.ps tmp/102nea1290537549.png",intern=TRUE))
character(0)
>
>
> proc.time()
user system elapsed
3.912 2.574 16.428