R version 3.0.2 (2013-09-25) -- "Frisbee Sailing"
Copyright (C) 2013 The R Foundation for Statistical Computing
Platform: i686-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(2.939217848
+ ,3
+ ,0.38
+ ,0.3
+ ,0.18
+ ,0.033717848
+ ,2.897889829
+ ,2
+ ,0.33
+ ,0.03
+ ,0.16
+ ,-0.049760171
+ ,2.821384807
+ ,3
+ ,0.14
+ ,0.47
+ ,0.92
+ ,0.000734807
+ ,2.976828188
+ ,1
+ ,0.13
+ ,0.79
+ ,0.11
+ ,0.028678188
+ ,2.98658494
+ ,1
+ ,0.24
+ ,0.81
+ ,0.43
+ ,0.04813494
+ ,3.08609928
+ ,1
+ ,0.96
+ ,0.42
+ ,0.09
+ ,0.02769928
+ ,2.970331644
+ ,2
+ ,0.4
+ ,0.17
+ ,0.34
+ ,0.030381644
+ ,3.001422707
+ ,3
+ ,0.98
+ ,0.75
+ ,0.69
+ ,0.070872707
+ ,2.896078452
+ ,2
+ ,0.01
+ ,0.72
+ ,0.25
+ ,0.010378452
+ ,3.048743903
+ ,1
+ ,0.57
+ ,0.77
+ ,0.65
+ ,0.085093903
+ ,2.951881247
+ ,1
+ ,0.08
+ ,0.06
+ ,0.13
+ ,-0.007118753
+ ,2.856323604
+ ,4
+ ,0.24
+ ,0.25
+ ,0.42
+ ,0.023173604
+ ,2.879911503
+ ,2
+ ,0.02
+ ,0.16
+ ,0.91
+ ,0.025211503
+ ,3.040257383
+ ,4
+ ,0.87
+ ,0.16
+ ,0.2
+ ,0.113857383
+ ,3.003077554
+ ,0
+ ,0.01
+ ,0.53
+ ,0.12
+ ,0.023527554
+ ,2.959791865
+ ,2
+ ,0.43
+ ,0.81
+ ,0.71
+ ,0.058141865
+ ,2.917990764
+ ,1
+ ,0.15
+ ,0.69
+ ,0.48
+ ,-0.009159236
+ ,2.818941548
+ ,2
+ ,0.65
+ ,0.86
+ ,0.85
+ ,-0.098058452
+ ,2.914778377
+ ,4
+ ,0.78
+ ,0.34
+ ,0.67
+ ,0.036578377
+ ,2.915101361
+ ,4
+ ,0.57
+ ,0.84
+ ,0.5
+ ,0.062701361
+ ,2.847462114
+ ,4
+ ,0.24
+ ,0.57
+ ,0.18
+ ,0.005512114
+ ,2.914117188
+ ,2
+ ,0.87
+ ,0.27
+ ,0.45
+ ,-0.072032812
+ ,2.923740777
+ ,4
+ ,0.5
+ ,0.58
+ ,0.66
+ ,0.084440777
+ ,2.913471302
+ ,4
+ ,0.26
+ ,0.04
+ ,0.31
+ ,0.064971302
+ ,2.768332582
+ ,4
+ ,0.13
+ ,0.78
+ ,0.32
+ ,-0.045367418
+ ,2.780356645
+ ,3
+ ,0.1
+ ,0.29
+ ,0.75
+ ,-0.051893355
+ ,2.872434637
+ ,2
+ ,0.16
+ ,0.43
+ ,0.69
+ ,-0.007715363
+ ,2.957005147
+ ,3
+ ,0.67
+ ,0.17
+ ,0.61
+ ,0.043555147
+ ,2.877572021
+ ,3
+ ,0.38
+ ,0.44
+ ,0.15
+ ,-0.026527979
+ ,2.892788186
+ ,0
+ ,0.03
+ ,0.14
+ ,0.83
+ ,-0.049211814
+ ,2.938104209
+ ,4
+ ,0.8
+ ,0.12
+ ,0.86
+ ,0.065304209
+ ,2.963092457
+ ,1
+ ,0.63
+ ,0.06
+ ,0.81
+ ,-0.014307543
+ ,2.947805711
+ ,2
+ ,0.31
+ ,0.1
+ ,0.18
+ ,0.005705711
+ ,2.929058982
+ ,2
+ ,0.1
+ ,0.33
+ ,0.78
+ ,0.059908982
+ ,2.798268951
+ ,3
+ ,0.22
+ ,0.76
+ ,0.81
+ ,-0.032431049
+ ,3.046956202
+ ,0
+ ,0.79
+ ,0.79
+ ,0.47
+ ,0.004806202
+ ,2.909470599
+ ,4
+ ,0.7
+ ,0.71
+ ,0.52
+ ,0.039620599
+ ,2.916937627
+ ,2
+ ,0.93
+ ,0.89
+ ,0.94
+ ,-0.026612373
+ ,2.936951397
+ ,0
+ ,0.19
+ ,0.35
+ ,0.67
+ ,-0.030198603
+ ,2.935183006
+ ,2
+ ,0.48
+ ,0.23
+ ,0.54
+ ,0.001133006
+ ,3.071599792
+ ,0
+ ,0.87
+ ,0.55
+ ,0.8
+ ,0.036949792
+ ,2.944966979
+ ,1
+ ,0.08
+ ,0.91
+ ,0.12
+ ,0.006516979
+ ,2.997397751
+ ,0
+ ,0.41
+ ,0.44
+ ,0.06
+ ,-0.036602249
+ ,3.048360762
+ ,1
+ ,0.8
+ ,0.5
+ ,0.18
+ ,0.017460762
+ ,2.859827106
+ ,1
+ ,0.13
+ ,0.64
+ ,0.37
+ ,-0.073872894
+ ,2.885496007
+ ,4
+ ,0.56
+ ,0.24
+ ,0.41
+ ,0.012996007
+ ,2.837868559
+ ,3
+ ,0.23
+ ,0.69
+ ,0.72
+ ,-0.002081441
+ ,2.879873437
+ ,2
+ ,0.34
+ ,0.83
+ ,0.78
+ ,-0.005576563
+ ,2.987063453
+ ,0
+ ,0.17
+ ,0.32
+ ,0.78
+ ,0.029263453
+ ,2.899181373
+ ,4
+ ,0.6
+ ,0.64
+ ,0.49
+ ,0.037481373
+ ,2.731576475
+ ,4
+ ,0.26
+ ,0.34
+ ,0.65
+ ,-0.085623525
+ ,2.908615135
+ ,2
+ ,0.9
+ ,0.17
+ ,0.48
+ ,-0.081534865
+ ,2.726225589
+ ,4
+ ,0.24
+ ,0.05
+ ,0.17
+ ,-0.129424411
+ ,2.883899312
+ ,3
+ ,0.91
+ ,0.77
+ ,0.81
+ ,-0.029350688
+ ,2.926219759
+ ,4
+ ,0.98
+ ,0.16
+ ,0.11
+ ,-0.019680241
+ ,2.933021763
+ ,1
+ ,0.54
+ ,0.48
+ ,0.74
+ ,-0.027978237
+ ,3.099391718
+ ,0
+ ,0.65
+ ,0.01
+ ,0.67
+ ,0.068541718
+ ,2.929512226
+ ,2
+ ,0.45
+ ,0.8
+ ,0.77
+ ,0.029412226
+ ,2.948690227
+ ,2
+ ,0.47
+ ,0.67
+ ,0.02
+ ,-0.009559773
+ ,2.997111915
+ ,4
+ ,0.57
+ ,0.85
+ ,0.31
+ ,0.131661915)
+ ,dim=c(6
+ ,60)
+ ,dimnames=list(c('echtscheiding'
+ ,'aantalkinderen'
+ ,'stresserendejob'
+ ,'seksleven'
+ ,'openkunnenpraten'
+ ,'storingsterm')
+ ,1:60))
> y <- array(NA,dim=c(6,60),dimnames=list(c('echtscheiding','aantalkinderen','stresserendejob','seksleven','openkunnenpraten','storingsterm'),1:60))
> for (i in 1:dim(x)[1])
+ {
+ for (j in 1:dim(x)[2])
+ {
+ y[i,j] <- as.numeric(x[i,j])
+ }
+ }
> par3 = 'No Linear Trend'
> par2 = 'Do not include Seasonal Dummies'
> par1 = '1'
> par3 <- 'No Linear Trend'
> par2 <- 'Do not include Seasonal Dummies'
> par1 <- '1'
> #'GNU S' R Code compiled by R2WASP v. 1.2.327 ()
> #Author: root
> #To cite this work: Wessa P., (2013), Multiple Regression (v1.0.29) in Free Statistics Software (v$_version), Office for Research Development and Education, URL http://www.wessa.net/rwasp_multipleregression.wasp/
> #Source of accompanying publication: Office for Research, Development, and Education
> #
> library(lattice)
> library(lmtest)
Loading required package: zoo
Attaching package: 'zoo'
The following objects are masked from 'package:base':
as.Date, 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
echtscheiding aantalkinderen stresserendejob seksleven openkunnenpraten
1 2.939218 3 0.38 0.30 0.18
2 2.897890 2 0.33 0.03 0.16
3 2.821385 3 0.14 0.47 0.92
4 2.976828 1 0.13 0.79 0.11
5 2.986585 1 0.24 0.81 0.43
6 3.086099 1 0.96 0.42 0.09
7 2.970332 2 0.40 0.17 0.34
8 3.001423 3 0.98 0.75 0.69
9 2.896078 2 0.01 0.72 0.25
10 3.048744 1 0.57 0.77 0.65
11 2.951881 1 0.08 0.06 0.13
12 2.856324 4 0.24 0.25 0.42
13 2.879912 2 0.02 0.16 0.91
14 3.040257 4 0.87 0.16 0.20
15 3.003078 0 0.01 0.53 0.12
16 2.959792 2 0.43 0.81 0.71
17 2.917991 1 0.15 0.69 0.48
18 2.818942 2 0.65 0.86 0.85
19 2.914778 4 0.78 0.34 0.67
20 2.915101 4 0.57 0.84 0.50
21 2.847462 4 0.24 0.57 0.18
22 2.914117 2 0.87 0.27 0.45
23 2.923741 4 0.50 0.58 0.66
24 2.913471 4 0.26 0.04 0.31
25 2.768333 4 0.13 0.78 0.32
26 2.780357 3 0.10 0.29 0.75
27 2.872435 2 0.16 0.43 0.69
28 2.957005 3 0.67 0.17 0.61
29 2.877572 3 0.38 0.44 0.15
30 2.892788 0 0.03 0.14 0.83
31 2.938104 4 0.80 0.12 0.86
32 2.963092 1 0.63 0.06 0.81
33 2.947806 2 0.31 0.10 0.18
34 2.929059 2 0.10 0.33 0.78
35 2.798269 3 0.22 0.76 0.81
36 3.046956 0 0.79 0.79 0.47
37 2.909471 4 0.70 0.71 0.52
38 2.916938 2 0.93 0.89 0.94
39 2.936951 0 0.19 0.35 0.67
40 2.935183 2 0.48 0.23 0.54
41 3.071600 0 0.87 0.55 0.80
42 2.944967 1 0.08 0.91 0.12
43 2.997398 0 0.41 0.44 0.06
44 3.048361 1 0.80 0.50 0.18
45 2.859827 1 0.13 0.64 0.37
46 2.885496 4 0.56 0.24 0.41
47 2.837869 3 0.23 0.69 0.72
48 2.879873 2 0.34 0.83 0.78
49 2.987063 0 0.17 0.32 0.78
50 2.899181 4 0.60 0.64 0.49
51 2.731576 4 0.26 0.34 0.65
52 2.908615 2 0.90 0.17 0.48
53 2.726226 4 0.24 0.05 0.17
54 2.883899 3 0.91 0.77 0.81
55 2.926220 4 0.98 0.16 0.11
56 2.933022 1 0.54 0.48 0.74
57 3.099392 0 0.65 0.01 0.67
58 2.929512 2 0.45 0.80 0.77
59 2.948690 2 0.47 0.67 0.02
60 2.997112 4 0.57 0.85 0.31
storingsterm
1 0.033717848
2 -0.049760171
3 0.000734807
4 0.028678188
5 0.048134940
6 0.027699280
7 0.030381644
8 0.070872707
9 0.010378452
10 0.085093903
11 -0.007118753
12 0.023173604
13 0.025211503
14 0.113857383
15 0.023527554
16 0.058141865
17 -0.009159236
18 -0.098058452
19 0.036578377
20 0.062701361
21 0.005512114
22 -0.072032812
23 0.084440777
24 0.064971302
25 -0.045367418
26 -0.051893355
27 -0.007715363
28 0.043555147
29 -0.026527979
30 -0.049211814
31 0.065304209
32 -0.014307543
33 0.005705711
34 0.059908982
35 -0.032431049
36 0.004806202
37 0.039620599
38 -0.026612373
39 -0.030198603
40 0.001133006
41 0.036949792
42 0.006516979
43 -0.036602249
44 0.017460762
45 -0.073872894
46 0.012996007
47 -0.002081441
48 -0.005576563
49 0.029263453
50 0.037481373
51 -0.085623525
52 -0.081534865
53 -0.129424411
54 -0.029350688
55 -0.019680241
56 -0.027978237
57 0.068541718
58 0.029412226
59 -0.009559773
60 0.131661915
> k <- length(x[1,])
> df <- as.data.frame(x)
> (mylm <- lm(df))
Call:
lm(formula = df)
Coefficients:
(Intercept) aantalkinderen stresserendejob seksleven
3.000 -0.040 0.120 -0.025
openkunnenpraten storingsterm
-0.070 1.000
> (mysum <- summary(mylm))
Call:
lm(formula = df)
Residuals:
Min 1Q Median 3Q Max
-2.975e-15 -7.611e-17 4.671e-17 1.339e-16 3.934e-16
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 3.000e+00 1.877e-16 1.598e+16 <2e-16 ***
aantalkinderen -4.000e-02 4.224e-17 -9.469e+14 <2e-16 ***
stresserendejob 1.200e-01 1.992e-16 6.025e+14 <2e-16 ***
seksleven -2.500e-02 2.068e-16 -1.209e+14 <2e-16 ***
openkunnenpraten -7.000e-02 2.107e-16 -3.322e+14 <2e-16 ***
storingsterm 1.000e+00 1.143e-15 8.747e+14 <2e-16 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 4.395e-16 on 54 degrees of freedom
Multiple R-squared: 1, Adjusted R-squared: 1
F-statistic: 3.91e+29 on 5 and 54 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,] 5.268505e-01 9.462990e-01 4.731495e-01
[2,] 7.141183e-01 5.717634e-01 2.858817e-01
[3,] 6.712416e-01 6.575168e-01 3.287584e-01
[4,] 9.992792e-01 1.441548e-03 7.207740e-04
[5,] 2.637436e-01 5.274873e-01 7.362564e-01
[6,] 1.000000e+00 2.688567e-13 1.344283e-13
[7,] 9.999002e-01 1.996528e-04 9.982640e-05
[8,] 4.506380e-01 9.012760e-01 5.493620e-01
[9,] 3.513520e-01 7.027040e-01 6.486480e-01
[10,] 9.999812e-01 3.764576e-05 1.882288e-05
[11,] 1.000000e+00 3.454259e-16 1.727130e-16
[12,] 4.436408e-02 8.872816e-02 9.556359e-01
[13,] 7.004082e-02 1.400816e-01 9.299592e-01
[14,] 1.163646e-01 2.327291e-01 8.836354e-01
[15,] 1.000000e+00 2.546338e-10 1.273169e-10
[16,] 8.739073e-01 2.521854e-01 1.260927e-01
[17,] 4.731432e-03 9.462865e-03 9.952686e-01
[18,] 9.127342e-01 1.745316e-01 8.726579e-02
[19,] 1.297954e-07 2.595907e-07 9.999999e-01
[20,] 1.397011e-01 2.794022e-01 8.602989e-01
[21,] 2.886476e-01 5.772952e-01 7.113524e-01
[22,] 8.252648e-01 3.494704e-01 1.747352e-01
[23,] 9.999991e-01 1.862428e-06 9.312142e-07
[24,] 9.999998e-01 3.101748e-07 1.550874e-07
[25,] 9.343539e-01 1.312921e-01 6.564607e-02
[26,] 7.196119e-01 5.607763e-01 2.803881e-01
[27,] 1.330167e-09 2.660334e-09 1.000000e+00
[28,] 9.785724e-01 4.285524e-02 2.142762e-02
[29,] 9.999958e-01 8.499362e-06 4.249681e-06
[30,] 4.517482e-04 9.034964e-04 9.995483e-01
[31,] 2.650364e-01 5.300728e-01 7.349636e-01
[32,] 1.974096e-01 3.948193e-01 8.025904e-01
[33,] 3.747778e-06 7.495556e-06 9.999963e-01
[34,] 3.675681e-02 7.351362e-02 9.632432e-01
[35,] 3.952751e-02 7.905502e-02 9.604725e-01
[36,] 3.438207e-01 6.876414e-01 6.561793e-01
[37,] 1.460234e-02 2.920469e-02 9.853977e-01
[38,] 5.872269e-08 1.174454e-07 9.999999e-01
[39,] 8.718467e-01 2.563066e-01 1.281533e-01
[40,] 5.322575e-01 9.354850e-01 4.677425e-01
[41,] 8.293204e-01 3.413591e-01 1.706796e-01
[42,] 2.760451e-01 5.520903e-01 7.239549e-01
[43,] 5.905730e-02 1.181146e-01 9.409427e-01
> postscript(file="/var/fisher/rcomp/tmp/1d0bl1384981363.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/fisher/rcomp/tmp/2oajs1384981363.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/fisher/rcomp/tmp/3e50b1384981363.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/fisher/rcomp/tmp/4thtd1384981363.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/fisher/rcomp/tmp/55al41384981363.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 = 60
Frequency = 1
1 2 3 4 5
-2.974938e-15 1.278601e-16 -2.018614e-16 -8.415659e-17 1.944755e-16
6 7 8 9 10
7.683692e-17 -6.733698e-17 8.092348e-17 -8.056644e-17 1.043863e-17
11 12 13 14 15
3.742454e-16 1.480405e-16 2.121837e-16 3.197276e-16 -7.521949e-17
16 17 18 19 20
8.296900e-17 -2.072315e-16 6.804128e-17 2.614575e-16 1.502917e-16
21 22 23 24 25
4.636217e-17 8.431899e-17 -1.308238e-17 2.765738e-16 3.934335e-16
26 27 28 29 30
1.728210e-18 2.181754e-17 -1.053593e-16 3.737818e-16 1.214324e-16
31 32 33 34 35
4.705927e-17 5.758857e-17 1.691288e-16 -6.071389e-17 -1.475968e-16
36 37 38 39 40
-8.498325e-17 1.291812e-16 -3.861568e-16 4.067968e-18 -2.891311e-17
41 42 43 44 45
2.114834e-17 2.345612e-16 -1.234287e-16 1.796291e-16 2.014808e-17
46 47 48 49 50
1.244522e-16 -7.878544e-17 1.972394e-17 1.027789e-16 -1.110378e-16
51 52 53 54 55
-7.901402e-17 -2.060430e-16 6.785949e-17 -2.883898e-17 6.046199e-17
56 57 58 59 60
8.118528e-17 3.819505e-17 -1.235903e-16 2.931870e-16 1.915578e-16
> postscript(file="/var/fisher/rcomp/tmp/6mudc1384981363.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 = 60
Frequency = 1
lag(myerror, k = 1) myerror
0 -2.974938e-15 NA
1 1.278601e-16 -2.974938e-15
2 -2.018614e-16 1.278601e-16
3 -8.415659e-17 -2.018614e-16
4 1.944755e-16 -8.415659e-17
5 7.683692e-17 1.944755e-16
6 -6.733698e-17 7.683692e-17
7 8.092348e-17 -6.733698e-17
8 -8.056644e-17 8.092348e-17
9 1.043863e-17 -8.056644e-17
10 3.742454e-16 1.043863e-17
11 1.480405e-16 3.742454e-16
12 2.121837e-16 1.480405e-16
13 3.197276e-16 2.121837e-16
14 -7.521949e-17 3.197276e-16
15 8.296900e-17 -7.521949e-17
16 -2.072315e-16 8.296900e-17
17 6.804128e-17 -2.072315e-16
18 2.614575e-16 6.804128e-17
19 1.502917e-16 2.614575e-16
20 4.636217e-17 1.502917e-16
21 8.431899e-17 4.636217e-17
22 -1.308238e-17 8.431899e-17
23 2.765738e-16 -1.308238e-17
24 3.934335e-16 2.765738e-16
25 1.728210e-18 3.934335e-16
26 2.181754e-17 1.728210e-18
27 -1.053593e-16 2.181754e-17
28 3.737818e-16 -1.053593e-16
29 1.214324e-16 3.737818e-16
30 4.705927e-17 1.214324e-16
31 5.758857e-17 4.705927e-17
32 1.691288e-16 5.758857e-17
33 -6.071389e-17 1.691288e-16
34 -1.475968e-16 -6.071389e-17
35 -8.498325e-17 -1.475968e-16
36 1.291812e-16 -8.498325e-17
37 -3.861568e-16 1.291812e-16
38 4.067968e-18 -3.861568e-16
39 -2.891311e-17 4.067968e-18
40 2.114834e-17 -2.891311e-17
41 2.345612e-16 2.114834e-17
42 -1.234287e-16 2.345612e-16
43 1.796291e-16 -1.234287e-16
44 2.014808e-17 1.796291e-16
45 1.244522e-16 2.014808e-17
46 -7.878544e-17 1.244522e-16
47 1.972394e-17 -7.878544e-17
48 1.027789e-16 1.972394e-17
49 -1.110378e-16 1.027789e-16
50 -7.901402e-17 -1.110378e-16
51 -2.060430e-16 -7.901402e-17
52 6.785949e-17 -2.060430e-16
53 -2.883898e-17 6.785949e-17
54 6.046199e-17 -2.883898e-17
55 8.118528e-17 6.046199e-17
56 3.819505e-17 8.118528e-17
57 -1.235903e-16 3.819505e-17
58 2.931870e-16 -1.235903e-16
59 1.915578e-16 2.931870e-16
60 NA 1.915578e-16
> dum1 <- dum[2:length(myerror),]
> dum1
lag(myerror, k = 1) myerror
[1,] 1.278601e-16 -2.974938e-15
[2,] -2.018614e-16 1.278601e-16
[3,] -8.415659e-17 -2.018614e-16
[4,] 1.944755e-16 -8.415659e-17
[5,] 7.683692e-17 1.944755e-16
[6,] -6.733698e-17 7.683692e-17
[7,] 8.092348e-17 -6.733698e-17
[8,] -8.056644e-17 8.092348e-17
[9,] 1.043863e-17 -8.056644e-17
[10,] 3.742454e-16 1.043863e-17
[11,] 1.480405e-16 3.742454e-16
[12,] 2.121837e-16 1.480405e-16
[13,] 3.197276e-16 2.121837e-16
[14,] -7.521949e-17 3.197276e-16
[15,] 8.296900e-17 -7.521949e-17
[16,] -2.072315e-16 8.296900e-17
[17,] 6.804128e-17 -2.072315e-16
[18,] 2.614575e-16 6.804128e-17
[19,] 1.502917e-16 2.614575e-16
[20,] 4.636217e-17 1.502917e-16
[21,] 8.431899e-17 4.636217e-17
[22,] -1.308238e-17 8.431899e-17
[23,] 2.765738e-16 -1.308238e-17
[24,] 3.934335e-16 2.765738e-16
[25,] 1.728210e-18 3.934335e-16
[26,] 2.181754e-17 1.728210e-18
[27,] -1.053593e-16 2.181754e-17
[28,] 3.737818e-16 -1.053593e-16
[29,] 1.214324e-16 3.737818e-16
[30,] 4.705927e-17 1.214324e-16
[31,] 5.758857e-17 4.705927e-17
[32,] 1.691288e-16 5.758857e-17
[33,] -6.071389e-17 1.691288e-16
[34,] -1.475968e-16 -6.071389e-17
[35,] -8.498325e-17 -1.475968e-16
[36,] 1.291812e-16 -8.498325e-17
[37,] -3.861568e-16 1.291812e-16
[38,] 4.067968e-18 -3.861568e-16
[39,] -2.891311e-17 4.067968e-18
[40,] 2.114834e-17 -2.891311e-17
[41,] 2.345612e-16 2.114834e-17
[42,] -1.234287e-16 2.345612e-16
[43,] 1.796291e-16 -1.234287e-16
[44,] 2.014808e-17 1.796291e-16
[45,] 1.244522e-16 2.014808e-17
[46,] -7.878544e-17 1.244522e-16
[47,] 1.972394e-17 -7.878544e-17
[48,] 1.027789e-16 1.972394e-17
[49,] -1.110378e-16 1.027789e-16
[50,] -7.901402e-17 -1.110378e-16
[51,] -2.060430e-16 -7.901402e-17
[52,] 6.785949e-17 -2.060430e-16
[53,] -2.883898e-17 6.785949e-17
[54,] 6.046199e-17 -2.883898e-17
[55,] 8.118528e-17 6.046199e-17
[56,] 3.819505e-17 8.118528e-17
[57,] -1.235903e-16 3.819505e-17
[58,] 2.931870e-16 -1.235903e-16
[59,] 1.915578e-16 2.931870e-16
> z <- as.data.frame(dum1)
> z
lag(myerror, k = 1) myerror
1 1.278601e-16 -2.974938e-15
2 -2.018614e-16 1.278601e-16
3 -8.415659e-17 -2.018614e-16
4 1.944755e-16 -8.415659e-17
5 7.683692e-17 1.944755e-16
6 -6.733698e-17 7.683692e-17
7 8.092348e-17 -6.733698e-17
8 -8.056644e-17 8.092348e-17
9 1.043863e-17 -8.056644e-17
10 3.742454e-16 1.043863e-17
11 1.480405e-16 3.742454e-16
12 2.121837e-16 1.480405e-16
13 3.197276e-16 2.121837e-16
14 -7.521949e-17 3.197276e-16
15 8.296900e-17 -7.521949e-17
16 -2.072315e-16 8.296900e-17
17 6.804128e-17 -2.072315e-16
18 2.614575e-16 6.804128e-17
19 1.502917e-16 2.614575e-16
20 4.636217e-17 1.502917e-16
21 8.431899e-17 4.636217e-17
22 -1.308238e-17 8.431899e-17
23 2.765738e-16 -1.308238e-17
24 3.934335e-16 2.765738e-16
25 1.728210e-18 3.934335e-16
26 2.181754e-17 1.728210e-18
27 -1.053593e-16 2.181754e-17
28 3.737818e-16 -1.053593e-16
29 1.214324e-16 3.737818e-16
30 4.705927e-17 1.214324e-16
31 5.758857e-17 4.705927e-17
32 1.691288e-16 5.758857e-17
33 -6.071389e-17 1.691288e-16
34 -1.475968e-16 -6.071389e-17
35 -8.498325e-17 -1.475968e-16
36 1.291812e-16 -8.498325e-17
37 -3.861568e-16 1.291812e-16
38 4.067968e-18 -3.861568e-16
39 -2.891311e-17 4.067968e-18
40 2.114834e-17 -2.891311e-17
41 2.345612e-16 2.114834e-17
42 -1.234287e-16 2.345612e-16
43 1.796291e-16 -1.234287e-16
44 2.014808e-17 1.796291e-16
45 1.244522e-16 2.014808e-17
46 -7.878544e-17 1.244522e-16
47 1.972394e-17 -7.878544e-17
48 1.027789e-16 1.972394e-17
49 -1.110378e-16 1.027789e-16
50 -7.901402e-17 -1.110378e-16
51 -2.060430e-16 -7.901402e-17
52 6.785949e-17 -2.060430e-16
53 -2.883898e-17 6.785949e-17
54 6.046199e-17 -2.883898e-17
55 8.118528e-17 6.046199e-17
56 3.819505e-17 8.118528e-17
57 -1.235903e-16 3.819505e-17
58 2.931870e-16 -1.235903e-16
59 1.915578e-16 2.931870e-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/fisher/rcomp/tmp/7y24u1384981363.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/fisher/rcomp/tmp/82pr51384981363.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/fisher/rcomp/tmp/9i5nc1384981363.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/fisher/rcomp/tmp/105n881384981363.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/fisher/rcomp/createtable file can be downloaded at http://www.wessa.net/cretab
> load(file="/var/fisher/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, signif(mysum$coefficients[i,1],6), 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/fisher/rcomp/tmp/11r9hu1384981363.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,signif(mysum$coefficients[i,1],6))
+ a<-table.element(a, signif(mysum$coefficients[i,2],6))
+ a<-table.element(a, signif(mysum$coefficients[i,3],4))
+ a<-table.element(a, signif(mysum$coefficients[i,4],6))
+ a<-table.element(a, signif(mysum$coefficients[i,4]/2,6))
+ a<-table.row.end(a)
+ }
> a<-table.end(a)
> table.save(a,file="/var/fisher/rcomp/tmp/12t6td1384981363.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, signif(sqrt(mysum$r.squared),6))
> a<-table.row.end(a)
> a<-table.row.start(a)
> a<-table.element(a, 'R-squared',1,TRUE)
> a<-table.element(a, signif(mysum$r.squared,6))
> a<-table.row.end(a)
> a<-table.row.start(a)
> a<-table.element(a, 'Adjusted R-squared',1,TRUE)
> a<-table.element(a, signif(mysum$adj.r.squared,6))
> a<-table.row.end(a)
> a<-table.row.start(a)
> a<-table.element(a, 'F-TEST (value)',1,TRUE)
> a<-table.element(a, signif(mysum$fstatistic[1],6))
> a<-table.row.end(a)
> a<-table.row.start(a)
> a<-table.element(a, 'F-TEST (DF numerator)',1,TRUE)
> a<-table.element(a, signif(mysum$fstatistic[2],6))
> a<-table.row.end(a)
> a<-table.row.start(a)
> a<-table.element(a, 'F-TEST (DF denominator)',1,TRUE)
> a<-table.element(a, signif(mysum$fstatistic[3],6))
> a<-table.row.end(a)
> a<-table.row.start(a)
> a<-table.element(a, 'p-value',1,TRUE)
> a<-table.element(a, signif(1-pf(mysum$fstatistic[1],mysum$fstatistic[2],mysum$fstatistic[3]),6))
> 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, signif(mysum$sigma,6))
> a<-table.row.end(a)
> a<-table.row.start(a)
> a<-table.element(a, 'Sum Squared Residuals',1,TRUE)
> a<-table.element(a, signif(sum(myerror*myerror),6))
> a<-table.row.end(a)
> a<-table.end(a)
> table.save(a,file="/var/fisher/rcomp/tmp/13uwr41384981364.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,signif(x[i],6))
+ a<-table.element(a,signif(x[i]-mysum$resid[i],6))
+ a<-table.element(a,signif(mysum$resid[i],6))
+ a<-table.row.end(a)
+ }
> a<-table.end(a)
> table.save(a,file="/var/fisher/rcomp/tmp/14dwsy1384981364.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,signif(gqarr[mypoint-kp3+1,1],6))
+ a<-table.element(a,signif(gqarr[mypoint-kp3+1,2],6))
+ a<-table.element(a,signif(gqarr[mypoint-kp3+1,3],6))
+ a<-table.row.end(a)
+ }
+ a<-table.end(a)
+ table.save(a,file="/var/fisher/rcomp/tmp/15yxo81384981364.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,signif(numsignificant1,6))
+ a<-table.element(a,signif(numsignificant1/numgqtests,6))
+ 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,signif(numsignificant5,6))
+ a<-table.element(a,signif(numsignificant5/numgqtests,6))
+ 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,signif(numsignificant10,6))
+ a<-table.element(a,signif(numsignificant10/numgqtests,6))
+ 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/fisher/rcomp/tmp/169kaz1384981364.tab")
+ }
>
> try(system("convert tmp/1d0bl1384981363.ps tmp/1d0bl1384981363.png",intern=TRUE))
character(0)
> try(system("convert tmp/2oajs1384981363.ps tmp/2oajs1384981363.png",intern=TRUE))
character(0)
> try(system("convert tmp/3e50b1384981363.ps tmp/3e50b1384981363.png",intern=TRUE))
character(0)
> try(system("convert tmp/4thtd1384981363.ps tmp/4thtd1384981363.png",intern=TRUE))
character(0)
> try(system("convert tmp/55al41384981363.ps tmp/55al41384981363.png",intern=TRUE))
character(0)
> try(system("convert tmp/6mudc1384981363.ps tmp/6mudc1384981363.png",intern=TRUE))
character(0)
> try(system("convert tmp/7y24u1384981363.ps tmp/7y24u1384981363.png",intern=TRUE))
character(0)
> try(system("convert tmp/82pr51384981363.ps tmp/82pr51384981363.png",intern=TRUE))
character(0)
> try(system("convert tmp/9i5nc1384981363.ps tmp/9i5nc1384981363.png",intern=TRUE))
character(0)
> try(system("convert tmp/105n881384981363.ps tmp/105n881384981363.png",intern=TRUE))
character(0)
>
>
> proc.time()
user system elapsed
5.768 1.238 7.013