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(20,115,21,112,22,109,23,101,24,107,25,99,26,120,27,133,28,128,29,135,30,160,31,144,32,168,33,112,34,170,35,110,36,150,37,144,38,124,39,156,40,160,41,133,42,177,43,186,44,190,45,133,46,154,47,179,48,166,49,177,50,160,51,150,52,177,53,159,54,147,55,161,56,168,57,120,58,170,59,180,60,123,61,167,62,150,63,123,64,144,65,154,66,178,67,181,68,100,69,132,70,120,71,146,72,160,73,120,74,158,75,148,76,190,77,138,78,189,79,170,80,156,81,177,82,165,83,181,84,147,85,163,86,169,87,173,88,133,89,156,90,179,91,126,92,163,93,156,94,136,95,183,96,189,97,170,98,167,99,189,100,178),dim=c(2,81),dimnames=list(c('Leeftijd','Bloeddruk'),1:81))
> y <- array(NA,dim=c(2,81),dimnames=list(c('Leeftijd','Bloeddruk'),1:81))
> 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
> 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
Leeftijd Bloeddruk t
1 20 115 1
2 21 112 2
3 22 109 3
4 23 101 4
5 24 107 5
6 25 99 6
7 26 120 7
8 27 133 8
9 28 128 9
10 29 135 10
11 30 160 11
12 31 144 12
13 32 168 13
14 33 112 14
15 34 170 15
16 35 110 16
17 36 150 17
18 37 144 18
19 38 124 19
20 39 156 20
21 40 160 21
22 41 133 22
23 42 177 23
24 43 186 24
25 44 190 25
26 45 133 26
27 46 154 27
28 47 179 28
29 48 166 29
30 49 177 30
31 50 160 31
32 51 150 32
33 52 177 33
34 53 159 34
35 54 147 35
36 55 161 36
37 56 168 37
38 57 120 38
39 58 170 39
40 59 180 40
41 60 123 41
42 61 167 42
43 62 150 43
44 63 123 44
45 64 144 45
46 65 154 46
47 66 178 47
48 67 181 48
49 68 100 49
50 69 132 50
51 70 120 51
52 71 146 52
53 72 160 53
54 73 120 54
55 74 158 55
56 75 148 56
57 76 190 57
58 77 138 58
59 78 189 59
60 79 170 60
61 80 156 61
62 81 177 62
63 82 165 63
64 83 181 64
65 84 147 65
66 85 163 66
67 86 169 67
68 87 173 68
69 88 133 69
70 89 156 70
71 90 179 71
72 91 126 72
73 92 163 73
74 93 156 74
75 94 136 75
76 95 183 76
77 96 189 77
78 97 170 78
79 98 167 79
80 99 189 80
81 100 178 81
> k <- length(x[1,])
> df <- as.data.frame(x)
> (mylm <- lm(df))
Call:
lm(formula = df)
Coefficients:
(Intercept) Bloeddruk t
19 0 1
> (mysum <- summary(mylm))
Call:
lm(formula = df)
Residuals:
Min 1Q Median 3Q Max
-1.365e-14 -1.164e-15 -1.052e-16 1.694e-15 4.815e-15
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 1.900e+01 1.664e-15 1.142e+16 <2e-16 ***
Bloeddruk 0.000e+00 1.195e-17 0.000e+00 1
t 1.000e+00 1.270e-17 7.875e+16 <2e-16 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 2.379e-15 on 78 degrees of freedom
Multiple R-squared: 1, Adjusted R-squared: 1
F-statistic: 3.912e+33 on 2 and 78 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,] 6.312015e-03 1.262403e-02 9.936880e-01
[2,] 1.602779e-04 3.205558e-04 9.998397e-01
[3,] 3.344273e-01 6.688546e-01 6.655727e-01
[4,] 6.641262e-02 1.328252e-01 9.335874e-01
[5,] 1.702064e-03 3.404128e-03 9.982979e-01
[6,] 2.470931e-02 4.941862e-02 9.752907e-01
[7,] 1.081477e-01 2.162954e-01 8.918523e-01
[8,] 1.772384e-03 3.544768e-03 9.982276e-01
[9,] 2.072695e-04 4.145390e-04 9.997927e-01
[10,] 6.647736e-08 1.329547e-07 9.999999e-01
[11,] 2.762610e-17 5.525220e-17 1.000000e+00
[12,] 2.999091e-02 5.998182e-02 9.700091e-01
[13,] 9.150567e-09 1.830113e-08 1.000000e+00
[14,] 3.103291e-15 6.206582e-15 1.000000e+00
[15,] 3.099642e-12 6.199283e-12 1.000000e+00
[16,] 1.495254e-10 2.990507e-10 1.000000e+00
[17,] 1.822861e-06 3.645723e-06 9.999982e-01
[18,] 6.947335e-03 1.389467e-02 9.930527e-01
[19,] 9.943095e-01 1.138102e-02 5.690508e-03
[20,] 6.412067e-15 1.282413e-14 1.000000e+00
[21,] 1.329035e-14 2.658069e-14 1.000000e+00
[22,] 9.228137e-06 1.845627e-05 9.999908e-01
[23,] 3.381755e-15 6.763510e-15 1.000000e+00
[24,] 4.028417e-07 8.056834e-07 9.999996e-01
[25,] 1.000000e+00 4.914524e-08 2.457262e-08
[26,] 2.102940e-02 4.205880e-02 9.789706e-01
[27,] 7.586784e-01 4.826433e-01 2.413216e-01
[28,] 5.684412e-01 8.631177e-01 4.315588e-01
[29,] 1.240165e-10 2.480330e-10 1.000000e+00
[30,] 1.206656e-12 2.413312e-12 1.000000e+00
[31,] 4.245459e-01 8.490918e-01 5.754541e-01
[32,] 5.865602e-12 1.173120e-11 1.000000e+00
[33,] 1.890938e-03 3.781877e-03 9.981091e-01
[34,] 2.414652e-15 4.829304e-15 1.000000e+00
[35,] 3.488787e-04 6.977573e-04 9.996511e-01
[36,] 4.188461e-09 8.376922e-09 1.000000e+00
[37,] 1.000000e+00 2.043829e-16 1.021914e-16
[38,] 6.565961e-08 1.313192e-07 9.999999e-01
[39,] 1.000000e+00 2.954944e-14 1.477472e-14
[40,] 7.999330e-04 1.599866e-03 9.992001e-01
[41,] 1.257162e-02 2.514324e-02 9.874284e-01
[42,] 9.999990e-01 1.969401e-06 9.847006e-07
[43,] 1.747980e-13 3.495959e-13 1.000000e+00
[44,] 2.326767e-03 4.653535e-03 9.976732e-01
[45,] 4.306873e-01 8.613746e-01 5.693127e-01
[46,] 4.940603e-04 9.881206e-04 9.995059e-01
[47,] 6.867435e-01 6.265130e-01 3.132565e-01
[48,] 1.356841e-01 2.713682e-01 8.643159e-01
[49,] 9.999998e-01 3.931971e-07 1.965985e-07
[50,] 2.879016e-06 5.758033e-06 9.999971e-01
[51,] 3.648419e-30 7.296837e-30 1.000000e+00
[52,] 3.726880e-11 7.453761e-11 1.000000e+00
[53,] 9.589347e-01 8.213055e-02 4.106528e-02
[54,] 8.040764e-01 3.918471e-01 1.959236e-01
[55,] 3.495483e-02 6.990966e-02 9.650452e-01
[56,] 9.990655e-01 1.869049e-03 9.345244e-04
[57,] 1.240646e-25 2.481292e-25 1.000000e+00
[58,] 2.918575e-02 5.837150e-02 9.708143e-01
[59,] 1.763429e-01 3.526857e-01 8.236571e-01
[60,] 1.000000e+00 5.194850e-11 2.597425e-11
[61,] 6.522027e-01 6.955946e-01 3.477973e-01
[62,] 4.296077e-01 8.592154e-01 5.703923e-01
[63,] 9.999868e-01 2.635017e-05 1.317509e-05
[64,] 1.872912e-01 3.745824e-01 8.127088e-01
[65,] 3.389244e-02 6.778488e-02 9.661076e-01
[66,] 2.214941e-01 4.429882e-01 7.785059e-01
[67,] 9.984312e-01 3.137556e-03 1.568778e-03
[68,] 3.168761e-10 6.337522e-10 1.000000e+00
[69,] 9.696953e-01 6.060946e-02 3.030473e-02
[70,] 6.665618e-01 6.668764e-01 3.334382e-01
> postscript(file="/var/wessaorg/rcomp/tmp/1l1ix1321385640.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/2lfzd1321385640.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/3ao8s1321385640.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/41tpz1321385640.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/5i0ej1321385640.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 = 81
Frequency = 1
1 2 3 4 5
-1.364582e-14 -9.116249e-16 1.457379e-15 -1.192148e-15 1.664764e-15
6 7 8 9 10
1.750084e-15 2.107564e-15 1.432981e-15 1.565373e-15 9.267503e-16
11 12 13 14 15
4.814560e-15 2.218963e-15 -1.052055e-16 2.533781e-15 -1.565534e-15
16 17 18 19 20
2.226659e-15 2.183239e-15 3.209080e-15 1.940379e-15 2.771574e-15
21 22 23 24 25
-2.288659e-15 -1.043715e-15 -3.645900e-15 -3.472805e-15 2.186414e-16
26 27 28 29 30
-1.239072e-15 -2.213209e-15 -1.056722e-15 -1.829324e-15 1.061950e-15
31 32 33 34 35
-6.375957e-16 7.601756e-16 -1.995705e-15 2.755978e-16 -4.589068e-16
36 37 38 39 40
1.813641e-16 -4.439415e-17 2.495821e-16 -1.026806e-15 -7.939295e-16
41 42 43 44 45
-9.859452e-16 1.376214e-16 6.259954e-16 -4.908865e-16 -1.909330e-15
46 47 48 49 50
-1.265324e-15 -1.164349e-15 -1.365122e-15 -1.622763e-15 -1.138946e-15
51 52 53 54 55
-2.803262e-15 -7.664689e-16 1.694394e-15 2.890763e-15 3.336885e-15
56 57 58 59 60
1.657690e-15 2.100511e-15 1.920740e-15 2.052338e-15 1.616346e-15
61 62 63 64 65
1.739795e-15 2.444653e-15 2.077910e-15 2.402328e-15 1.437700e-15
66 67 68 69 70
1.762985e-15 1.067184e-15 1.207651e-15 -1.497374e-15 -2.926358e-15
71 72 73 74 75
-1.122095e-16 -1.151036e-15 -6.345912e-16 -6.217622e-16 -5.889875e-16
76 77 78 79 80
-9.579161e-16 -7.013548e-16 -1.563764e-15 -1.317013e-15 -2.062961e-16
81
-2.765796e-15
> postscript(file="/var/wessaorg/rcomp/tmp/655nj1321385640.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 = 81
Frequency = 1
lag(myerror, k = 1) myerror
0 -1.364582e-14 NA
1 -9.116249e-16 -1.364582e-14
2 1.457379e-15 -9.116249e-16
3 -1.192148e-15 1.457379e-15
4 1.664764e-15 -1.192148e-15
5 1.750084e-15 1.664764e-15
6 2.107564e-15 1.750084e-15
7 1.432981e-15 2.107564e-15
8 1.565373e-15 1.432981e-15
9 9.267503e-16 1.565373e-15
10 4.814560e-15 9.267503e-16
11 2.218963e-15 4.814560e-15
12 -1.052055e-16 2.218963e-15
13 2.533781e-15 -1.052055e-16
14 -1.565534e-15 2.533781e-15
15 2.226659e-15 -1.565534e-15
16 2.183239e-15 2.226659e-15
17 3.209080e-15 2.183239e-15
18 1.940379e-15 3.209080e-15
19 2.771574e-15 1.940379e-15
20 -2.288659e-15 2.771574e-15
21 -1.043715e-15 -2.288659e-15
22 -3.645900e-15 -1.043715e-15
23 -3.472805e-15 -3.645900e-15
24 2.186414e-16 -3.472805e-15
25 -1.239072e-15 2.186414e-16
26 -2.213209e-15 -1.239072e-15
27 -1.056722e-15 -2.213209e-15
28 -1.829324e-15 -1.056722e-15
29 1.061950e-15 -1.829324e-15
30 -6.375957e-16 1.061950e-15
31 7.601756e-16 -6.375957e-16
32 -1.995705e-15 7.601756e-16
33 2.755978e-16 -1.995705e-15
34 -4.589068e-16 2.755978e-16
35 1.813641e-16 -4.589068e-16
36 -4.439415e-17 1.813641e-16
37 2.495821e-16 -4.439415e-17
38 -1.026806e-15 2.495821e-16
39 -7.939295e-16 -1.026806e-15
40 -9.859452e-16 -7.939295e-16
41 1.376214e-16 -9.859452e-16
42 6.259954e-16 1.376214e-16
43 -4.908865e-16 6.259954e-16
44 -1.909330e-15 -4.908865e-16
45 -1.265324e-15 -1.909330e-15
46 -1.164349e-15 -1.265324e-15
47 -1.365122e-15 -1.164349e-15
48 -1.622763e-15 -1.365122e-15
49 -1.138946e-15 -1.622763e-15
50 -2.803262e-15 -1.138946e-15
51 -7.664689e-16 -2.803262e-15
52 1.694394e-15 -7.664689e-16
53 2.890763e-15 1.694394e-15
54 3.336885e-15 2.890763e-15
55 1.657690e-15 3.336885e-15
56 2.100511e-15 1.657690e-15
57 1.920740e-15 2.100511e-15
58 2.052338e-15 1.920740e-15
59 1.616346e-15 2.052338e-15
60 1.739795e-15 1.616346e-15
61 2.444653e-15 1.739795e-15
62 2.077910e-15 2.444653e-15
63 2.402328e-15 2.077910e-15
64 1.437700e-15 2.402328e-15
65 1.762985e-15 1.437700e-15
66 1.067184e-15 1.762985e-15
67 1.207651e-15 1.067184e-15
68 -1.497374e-15 1.207651e-15
69 -2.926358e-15 -1.497374e-15
70 -1.122095e-16 -2.926358e-15
71 -1.151036e-15 -1.122095e-16
72 -6.345912e-16 -1.151036e-15
73 -6.217622e-16 -6.345912e-16
74 -5.889875e-16 -6.217622e-16
75 -9.579161e-16 -5.889875e-16
76 -7.013548e-16 -9.579161e-16
77 -1.563764e-15 -7.013548e-16
78 -1.317013e-15 -1.563764e-15
79 -2.062961e-16 -1.317013e-15
80 -2.765796e-15 -2.062961e-16
81 NA -2.765796e-15
> dum1 <- dum[2:length(myerror),]
> dum1
lag(myerror, k = 1) myerror
[1,] -9.116249e-16 -1.364582e-14
[2,] 1.457379e-15 -9.116249e-16
[3,] -1.192148e-15 1.457379e-15
[4,] 1.664764e-15 -1.192148e-15
[5,] 1.750084e-15 1.664764e-15
[6,] 2.107564e-15 1.750084e-15
[7,] 1.432981e-15 2.107564e-15
[8,] 1.565373e-15 1.432981e-15
[9,] 9.267503e-16 1.565373e-15
[10,] 4.814560e-15 9.267503e-16
[11,] 2.218963e-15 4.814560e-15
[12,] -1.052055e-16 2.218963e-15
[13,] 2.533781e-15 -1.052055e-16
[14,] -1.565534e-15 2.533781e-15
[15,] 2.226659e-15 -1.565534e-15
[16,] 2.183239e-15 2.226659e-15
[17,] 3.209080e-15 2.183239e-15
[18,] 1.940379e-15 3.209080e-15
[19,] 2.771574e-15 1.940379e-15
[20,] -2.288659e-15 2.771574e-15
[21,] -1.043715e-15 -2.288659e-15
[22,] -3.645900e-15 -1.043715e-15
[23,] -3.472805e-15 -3.645900e-15
[24,] 2.186414e-16 -3.472805e-15
[25,] -1.239072e-15 2.186414e-16
[26,] -2.213209e-15 -1.239072e-15
[27,] -1.056722e-15 -2.213209e-15
[28,] -1.829324e-15 -1.056722e-15
[29,] 1.061950e-15 -1.829324e-15
[30,] -6.375957e-16 1.061950e-15
[31,] 7.601756e-16 -6.375957e-16
[32,] -1.995705e-15 7.601756e-16
[33,] 2.755978e-16 -1.995705e-15
[34,] -4.589068e-16 2.755978e-16
[35,] 1.813641e-16 -4.589068e-16
[36,] -4.439415e-17 1.813641e-16
[37,] 2.495821e-16 -4.439415e-17
[38,] -1.026806e-15 2.495821e-16
[39,] -7.939295e-16 -1.026806e-15
[40,] -9.859452e-16 -7.939295e-16
[41,] 1.376214e-16 -9.859452e-16
[42,] 6.259954e-16 1.376214e-16
[43,] -4.908865e-16 6.259954e-16
[44,] -1.909330e-15 -4.908865e-16
[45,] -1.265324e-15 -1.909330e-15
[46,] -1.164349e-15 -1.265324e-15
[47,] -1.365122e-15 -1.164349e-15
[48,] -1.622763e-15 -1.365122e-15
[49,] -1.138946e-15 -1.622763e-15
[50,] -2.803262e-15 -1.138946e-15
[51,] -7.664689e-16 -2.803262e-15
[52,] 1.694394e-15 -7.664689e-16
[53,] 2.890763e-15 1.694394e-15
[54,] 3.336885e-15 2.890763e-15
[55,] 1.657690e-15 3.336885e-15
[56,] 2.100511e-15 1.657690e-15
[57,] 1.920740e-15 2.100511e-15
[58,] 2.052338e-15 1.920740e-15
[59,] 1.616346e-15 2.052338e-15
[60,] 1.739795e-15 1.616346e-15
[61,] 2.444653e-15 1.739795e-15
[62,] 2.077910e-15 2.444653e-15
[63,] 2.402328e-15 2.077910e-15
[64,] 1.437700e-15 2.402328e-15
[65,] 1.762985e-15 1.437700e-15
[66,] 1.067184e-15 1.762985e-15
[67,] 1.207651e-15 1.067184e-15
[68,] -1.497374e-15 1.207651e-15
[69,] -2.926358e-15 -1.497374e-15
[70,] -1.122095e-16 -2.926358e-15
[71,] -1.151036e-15 -1.122095e-16
[72,] -6.345912e-16 -1.151036e-15
[73,] -6.217622e-16 -6.345912e-16
[74,] -5.889875e-16 -6.217622e-16
[75,] -9.579161e-16 -5.889875e-16
[76,] -7.013548e-16 -9.579161e-16
[77,] -1.563764e-15 -7.013548e-16
[78,] -1.317013e-15 -1.563764e-15
[79,] -2.062961e-16 -1.317013e-15
[80,] -2.765796e-15 -2.062961e-16
> z <- as.data.frame(dum1)
> z
lag(myerror, k = 1) myerror
1 -9.116249e-16 -1.364582e-14
2 1.457379e-15 -9.116249e-16
3 -1.192148e-15 1.457379e-15
4 1.664764e-15 -1.192148e-15
5 1.750084e-15 1.664764e-15
6 2.107564e-15 1.750084e-15
7 1.432981e-15 2.107564e-15
8 1.565373e-15 1.432981e-15
9 9.267503e-16 1.565373e-15
10 4.814560e-15 9.267503e-16
11 2.218963e-15 4.814560e-15
12 -1.052055e-16 2.218963e-15
13 2.533781e-15 -1.052055e-16
14 -1.565534e-15 2.533781e-15
15 2.226659e-15 -1.565534e-15
16 2.183239e-15 2.226659e-15
17 3.209080e-15 2.183239e-15
18 1.940379e-15 3.209080e-15
19 2.771574e-15 1.940379e-15
20 -2.288659e-15 2.771574e-15
21 -1.043715e-15 -2.288659e-15
22 -3.645900e-15 -1.043715e-15
23 -3.472805e-15 -3.645900e-15
24 2.186414e-16 -3.472805e-15
25 -1.239072e-15 2.186414e-16
26 -2.213209e-15 -1.239072e-15
27 -1.056722e-15 -2.213209e-15
28 -1.829324e-15 -1.056722e-15
29 1.061950e-15 -1.829324e-15
30 -6.375957e-16 1.061950e-15
31 7.601756e-16 -6.375957e-16
32 -1.995705e-15 7.601756e-16
33 2.755978e-16 -1.995705e-15
34 -4.589068e-16 2.755978e-16
35 1.813641e-16 -4.589068e-16
36 -4.439415e-17 1.813641e-16
37 2.495821e-16 -4.439415e-17
38 -1.026806e-15 2.495821e-16
39 -7.939295e-16 -1.026806e-15
40 -9.859452e-16 -7.939295e-16
41 1.376214e-16 -9.859452e-16
42 6.259954e-16 1.376214e-16
43 -4.908865e-16 6.259954e-16
44 -1.909330e-15 -4.908865e-16
45 -1.265324e-15 -1.909330e-15
46 -1.164349e-15 -1.265324e-15
47 -1.365122e-15 -1.164349e-15
48 -1.622763e-15 -1.365122e-15
49 -1.138946e-15 -1.622763e-15
50 -2.803262e-15 -1.138946e-15
51 -7.664689e-16 -2.803262e-15
52 1.694394e-15 -7.664689e-16
53 2.890763e-15 1.694394e-15
54 3.336885e-15 2.890763e-15
55 1.657690e-15 3.336885e-15
56 2.100511e-15 1.657690e-15
57 1.920740e-15 2.100511e-15
58 2.052338e-15 1.920740e-15
59 1.616346e-15 2.052338e-15
60 1.739795e-15 1.616346e-15
61 2.444653e-15 1.739795e-15
62 2.077910e-15 2.444653e-15
63 2.402328e-15 2.077910e-15
64 1.437700e-15 2.402328e-15
65 1.762985e-15 1.437700e-15
66 1.067184e-15 1.762985e-15
67 1.207651e-15 1.067184e-15
68 -1.497374e-15 1.207651e-15
69 -2.926358e-15 -1.497374e-15
70 -1.122095e-16 -2.926358e-15
71 -1.151036e-15 -1.122095e-16
72 -6.345912e-16 -1.151036e-15
73 -6.217622e-16 -6.345912e-16
74 -5.889875e-16 -6.217622e-16
75 -9.579161e-16 -5.889875e-16
76 -7.013548e-16 -9.579161e-16
77 -1.563764e-15 -7.013548e-16
78 -1.317013e-15 -1.563764e-15
79 -2.062961e-16 -1.317013e-15
80 -2.765796e-15 -2.062961e-16
> 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/7s1gs1321385640.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/8jen41321385640.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/9f5201321385641.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/10l1ev1321385641.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/11pmor1321385641.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/12gery1321385641.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/13e5ik1321385641.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/14b6441321385641.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/157kwq1321385641.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/16yjvd1321385641.tab")
+ }
>
> try(system("convert tmp/1l1ix1321385640.ps tmp/1l1ix1321385640.png",intern=TRUE))
character(0)
> try(system("convert tmp/2lfzd1321385640.ps tmp/2lfzd1321385640.png",intern=TRUE))
character(0)
> try(system("convert tmp/3ao8s1321385640.ps tmp/3ao8s1321385640.png",intern=TRUE))
character(0)
> try(system("convert tmp/41tpz1321385640.ps tmp/41tpz1321385640.png",intern=TRUE))
character(0)
> try(system("convert tmp/5i0ej1321385640.ps tmp/5i0ej1321385640.png",intern=TRUE))
character(0)
> try(system("convert tmp/655nj1321385640.ps tmp/655nj1321385640.png",intern=TRUE))
character(0)
> try(system("convert tmp/7s1gs1321385640.ps tmp/7s1gs1321385640.png",intern=TRUE))
character(0)
> try(system("convert tmp/8jen41321385640.ps tmp/8jen41321385640.png",intern=TRUE))
character(0)
> try(system("convert tmp/9f5201321385641.ps tmp/9f5201321385641.png",intern=TRUE))
character(0)
> try(system("convert tmp/10l1ev1321385641.ps tmp/10l1ev1321385641.png",intern=TRUE))
character(0)
>
>
> proc.time()
user system elapsed
3.557 0.524 4.119