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(13328,0,0,12873,0,0,14000,0,0,13477,0,0,14237,0,0,13674,0,0,13529,0,0,14058,0,0,12975,0,0,14326,0,0,14008,0,0,16193,1,0,14483,0,0,14011,0,0,15057,0,0,14884,0,0,15414,0,0,14440,0,0,14900,0,0,15074,0,0,14442,0,0,15307,0,0,14938,0,0,17193,0,0,15528,0,0,14765,0,0,15838,0,0,15723,0,0,16150,0,0,15486,0,0,15986,0,1,15983,0,1,15692,0,1,16490,0,1,15686,0,1,18897,0,1,16316,0,1,15636,0,1,17163,0,1,16534,0,1,16518,0,1,16375,0,1,16290,0,1,16352,0,1,15943,0,1,16362,0,1,16393,0,1,19051,0,1,16747,0,1,16320,0,1,17910,0,1,16961,0,1,17480,0,1,17049,0,1,16879,0,1,17473,0,1,16998,0,1,17307,0,1,17418,0,1,20169,0,1,17871,0,1,17226,0,1,19062,0,1,17804,0,1,19100,0,1,18522,0,1,18060,0,1,18869,0,1,18127,0,1,18871,0,1,18890,0,1,21263,0,1,19547,0,1,18450,0,1,20254,0,1,19240,0,1,20216,0,1,19420,0,1,19415,0,1,20018,0,1,18652,0,1,19978,0,1,19509,0,1,21971,0,1),dim=c(3,84),dimnames=list(c('HPC','Dummies1','Dummies2'),1:84))
> y <- array(NA,dim=c(3,84),dimnames=list(c('HPC','Dummies1','Dummies2'),1:84))
> for (i in 1:dim(x)[1])
+ {
+ for (j in 1:dim(x)[2])
+ {
+ y[i,j] <- as.numeric(x[i,j])
+ }
+ }
> par3 = 'Linear Trend'
> par2 = 'Include Monthly Dummies'
> par1 = '1'
> par3 <- 'Linear Trend'
> par2 <- 'Include Monthly 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
HPC Dummies1 Dummies2 M1 M2 M3 M4 M5 M6 M7 M8 M9 M10 M11 t
1 13328 0 0 1 0 0 0 0 0 0 0 0 0 0 1
2 12873 0 0 0 1 0 0 0 0 0 0 0 0 0 2
3 14000 0 0 0 0 1 0 0 0 0 0 0 0 0 3
4 13477 0 0 0 0 0 1 0 0 0 0 0 0 0 4
5 14237 0 0 0 0 0 0 1 0 0 0 0 0 0 5
6 13674 0 0 0 0 0 0 0 1 0 0 0 0 0 6
7 13529 0 0 0 0 0 0 0 0 1 0 0 0 0 7
8 14058 0 0 0 0 0 0 0 0 0 1 0 0 0 8
9 12975 0 0 0 0 0 0 0 0 0 0 1 0 0 9
10 14326 0 0 0 0 0 0 0 0 0 0 0 1 0 10
11 14008 0 0 0 0 0 0 0 0 0 0 0 0 1 11
12 16193 1 0 0 0 0 0 0 0 0 0 0 0 0 12
13 14483 0 0 1 0 0 0 0 0 0 0 0 0 0 13
14 14011 0 0 0 1 0 0 0 0 0 0 0 0 0 14
15 15057 0 0 0 0 1 0 0 0 0 0 0 0 0 15
16 14884 0 0 0 0 0 1 0 0 0 0 0 0 0 16
17 15414 0 0 0 0 0 0 1 0 0 0 0 0 0 17
18 14440 0 0 0 0 0 0 0 1 0 0 0 0 0 18
19 14900 0 0 0 0 0 0 0 0 1 0 0 0 0 19
20 15074 0 0 0 0 0 0 0 0 0 1 0 0 0 20
21 14442 0 0 0 0 0 0 0 0 0 0 1 0 0 21
22 15307 0 0 0 0 0 0 0 0 0 0 0 1 0 22
23 14938 0 0 0 0 0 0 0 0 0 0 0 0 1 23
24 17193 0 0 0 0 0 0 0 0 0 0 0 0 0 24
25 15528 0 0 1 0 0 0 0 0 0 0 0 0 0 25
26 14765 0 0 0 1 0 0 0 0 0 0 0 0 0 26
27 15838 0 0 0 0 1 0 0 0 0 0 0 0 0 27
28 15723 0 0 0 0 0 1 0 0 0 0 0 0 0 28
29 16150 0 0 0 0 0 0 1 0 0 0 0 0 0 29
30 15486 0 0 0 0 0 0 0 1 0 0 0 0 0 30
31 15986 0 1 0 0 0 0 0 0 1 0 0 0 0 31
32 15983 0 1 0 0 0 0 0 0 0 1 0 0 0 32
33 15692 0 1 0 0 0 0 0 0 0 0 1 0 0 33
34 16490 0 1 0 0 0 0 0 0 0 0 0 1 0 34
35 15686 0 1 0 0 0 0 0 0 0 0 0 0 1 35
36 18897 0 1 0 0 0 0 0 0 0 0 0 0 0 36
37 16316 0 1 1 0 0 0 0 0 0 0 0 0 0 37
38 15636 0 1 0 1 0 0 0 0 0 0 0 0 0 38
39 17163 0 1 0 0 1 0 0 0 0 0 0 0 0 39
40 16534 0 1 0 0 0 1 0 0 0 0 0 0 0 40
41 16518 0 1 0 0 0 0 1 0 0 0 0 0 0 41
42 16375 0 1 0 0 0 0 0 1 0 0 0 0 0 42
43 16290 0 1 0 0 0 0 0 0 1 0 0 0 0 43
44 16352 0 1 0 0 0 0 0 0 0 1 0 0 0 44
45 15943 0 1 0 0 0 0 0 0 0 0 1 0 0 45
46 16362 0 1 0 0 0 0 0 0 0 0 0 1 0 46
47 16393 0 1 0 0 0 0 0 0 0 0 0 0 1 47
48 19051 0 1 0 0 0 0 0 0 0 0 0 0 0 48
49 16747 0 1 1 0 0 0 0 0 0 0 0 0 0 49
50 16320 0 1 0 1 0 0 0 0 0 0 0 0 0 50
51 17910 0 1 0 0 1 0 0 0 0 0 0 0 0 51
52 16961 0 1 0 0 0 1 0 0 0 0 0 0 0 52
53 17480 0 1 0 0 0 0 1 0 0 0 0 0 0 53
54 17049 0 1 0 0 0 0 0 1 0 0 0 0 0 54
55 16879 0 1 0 0 0 0 0 0 1 0 0 0 0 55
56 17473 0 1 0 0 0 0 0 0 0 1 0 0 0 56
57 16998 0 1 0 0 0 0 0 0 0 0 1 0 0 57
58 17307 0 1 0 0 0 0 0 0 0 0 0 1 0 58
59 17418 0 1 0 0 0 0 0 0 0 0 0 0 1 59
60 20169 0 1 0 0 0 0 0 0 0 0 0 0 0 60
61 17871 0 1 1 0 0 0 0 0 0 0 0 0 0 61
62 17226 0 1 0 1 0 0 0 0 0 0 0 0 0 62
63 19062 0 1 0 0 1 0 0 0 0 0 0 0 0 63
64 17804 0 1 0 0 0 1 0 0 0 0 0 0 0 64
65 19100 0 1 0 0 0 0 1 0 0 0 0 0 0 65
66 18522 0 1 0 0 0 0 0 1 0 0 0 0 0 66
67 18060 0 1 0 0 0 0 0 0 1 0 0 0 0 67
68 18869 0 1 0 0 0 0 0 0 0 1 0 0 0 68
69 18127 0 1 0 0 0 0 0 0 0 0 1 0 0 69
70 18871 0 1 0 0 0 0 0 0 0 0 0 1 0 70
71 18890 0 1 0 0 0 0 0 0 0 0 0 0 1 71
72 21263 0 1 0 0 0 0 0 0 0 0 0 0 0 72
73 19547 0 1 1 0 0 0 0 0 0 0 0 0 0 73
74 18450 0 1 0 1 0 0 0 0 0 0 0 0 0 74
75 20254 0 1 0 0 1 0 0 0 0 0 0 0 0 75
76 19240 0 1 0 0 0 1 0 0 0 0 0 0 0 76
77 20216 0 1 0 0 0 0 1 0 0 0 0 0 0 77
78 19420 0 1 0 0 0 0 0 1 0 0 0 0 0 78
79 19415 0 1 0 0 0 0 0 0 1 0 0 0 0 79
80 20018 0 1 0 0 0 0 0 0 0 1 0 0 0 80
81 18652 0 1 0 0 0 0 0 0 0 0 1 0 0 81
82 19978 0 1 0 0 0 0 0 0 0 0 0 1 0 82
83 19509 0 1 0 0 0 0 0 0 0 0 0 0 1 83
84 21971 0 1 0 0 0 0 0 0 0 0 0 0 0 84
> k <- length(x[1,])
> df <- as.data.frame(x)
> (mylm <- lm(df))
Call:
lm(formula = df)
Coefficients:
(Intercept) Dummies1 Dummies2 M1 M2 M3
15561.88 -338.96 -203.85 -2176.45 -2905.72 -1557.56
M4 M5 M6 M7 M8 M9
-2304.26 -1743.38 -2416.94 -2455.37 -2140.78 -2935.62
M10 M11 t
-2186.17 -2524.01 80.84
> (mysum <- summary(mylm))
Call:
lm(formula = df)
Residuals:
Min 1Q Median 3Q Max
-553.54 -185.99 19.09 169.84 628.75
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 15561.884 140.078 111.094 < 2e-16 ***
Dummies1 -338.956 321.595 -1.054 0.2956
Dummies2 -203.848 120.720 -1.689 0.0958 .
M1 -2176.454 164.400 -13.239 < 2e-16 ***
M2 -2905.722 164.208 -17.695 < 2e-16 ***
M3 -1557.561 164.049 -9.494 3.79e-14 ***
M4 -2304.257 163.926 -14.057 < 2e-16 ***
M5 -1743.382 163.837 -10.641 3.38e-16 ***
M6 -2416.936 163.782 -14.757 < 2e-16 ***
M7 -2455.369 163.614 -15.007 < 2e-16 ***
M8 -2140.779 163.422 -13.100 < 2e-16 ***
M9 -2935.619 163.265 -17.981 < 2e-16 ***
M10 -2186.172 163.142 -13.400 < 2e-16 ***
M11 -2524.011 163.053 -15.480 < 2e-16 ***
t 80.839 2.383 33.926 < 2e-16 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 292.6 on 69 degrees of freedom
Multiple R-squared: 0.9836, Adjusted R-squared: 0.9803
F-statistic: 295.3 on 14 and 69 DF, p-value: < 2.2e-16
> if (n > n25) {
+ kp3 <- k + 3
+ nmkm3 <- n - k - 3
+ gqarr <- array(NA, dim=c(nmkm3-kp3+1,3))
+ numgqtests <- 0
+ numsignificant1 <- 0
+ numsignificant5 <- 0
+ numsignificant10 <- 0
+ for (mypoint in kp3:nmkm3) {
+ j <- 0
+ numgqtests <- numgqtests + 1
+ for (myalt in c('greater', 'two.sided', 'less')) {
+ j <- j + 1
+ gqarr[mypoint-kp3+1,j] <- gqtest(mylm, point=mypoint, alternative=myalt)$p.value
+ }
+ if (gqarr[mypoint-kp3+1,2] < 0.01) numsignificant1 <- numsignificant1 + 1
+ if (gqarr[mypoint-kp3+1,2] < 0.05) numsignificant5 <- numsignificant5 + 1
+ if (gqarr[mypoint-kp3+1,2] < 0.10) numsignificant10 <- numsignificant10 + 1
+ }
+ gqarr
+ }
[,1] [,2] [,3]
[1,] 0.211260621 0.422521243 0.7887394
[2,] 0.142104874 0.284209749 0.8578951
[3,] 0.074134959 0.148269917 0.9258650
[4,] 0.065270334 0.130540668 0.9347297
[5,] 0.039145994 0.078291989 0.9608540
[6,] 0.024802375 0.049604751 0.9751976
[7,] 0.011881464 0.023762929 0.9881185
[8,] 0.005463828 0.010927656 0.9945362
[9,] 0.007189797 0.014379594 0.9928102
[10,] 0.008275011 0.016550022 0.9917250
[11,] 0.004168471 0.008336941 0.9958315
[12,] 0.003103776 0.006207553 0.9968962
[13,] 0.001553304 0.003106609 0.9984467
[14,] 0.001258451 0.002516902 0.9987415
[15,] 0.001122525 0.002245050 0.9988775
[16,] 0.002590701 0.005181401 0.9974093
[17,] 0.004420932 0.008841863 0.9955791
[18,] 0.010101349 0.020202698 0.9898987
[19,] 0.070669242 0.141338484 0.9293308
[20,] 0.086952127 0.173904253 0.9130479
[21,] 0.131716419 0.263432838 0.8682836
[22,] 0.110388633 0.220777266 0.8896114
[23,] 0.239660363 0.479320726 0.7603396
[24,] 0.527143937 0.945712126 0.4728561
[25,] 0.500761020 0.998477960 0.4992390
[26,] 0.633142407 0.733715186 0.3668576
[27,] 0.714655861 0.570688279 0.2853441
[28,] 0.725795106 0.548409787 0.2742049
[29,] 0.846956243 0.306087514 0.1530438
[30,] 0.815152879 0.369694242 0.1848471
[31,] 0.801888140 0.396223720 0.1981119
[32,] 0.792845747 0.414308505 0.2071543
[33,] 0.759351651 0.481296698 0.2406483
[34,] 0.698831470 0.602337060 0.3011685
[35,] 0.675787216 0.648425568 0.3242128
[36,] 0.711344076 0.577311848 0.2886559
[37,] 0.658195595 0.683608809 0.3418044
[38,] 0.625687668 0.748624665 0.3743123
[39,] 0.591417555 0.817164891 0.4085824
[40,] 0.554002858 0.891994285 0.4459971
[41,] 0.639225614 0.721548771 0.3607744
[42,] 0.564964489 0.870071023 0.4350355
[43,] 0.482560711 0.965121422 0.5174393
[44,] 0.624200705 0.751598590 0.3757993
[45,] 0.539577816 0.920844367 0.4604222
[46,] 0.486097925 0.972195849 0.5139021
[47,] 0.545321473 0.909357053 0.4546785
[48,] 0.468874263 0.937748526 0.5311257
[49,] 0.333801853 0.667603706 0.6661981
> postscript(file="/var/wessaorg/rcomp/tmp/1vhk41385408508.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/2pyn41385408508.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/3haag1385408508.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/4td8r1385408508.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/56gy61385408508.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 = 84
Frequency = 1
1 2 3 4 5
-1.382694e+02 5.515918e+01 -2.468408e+02 -1.039837e+02 1.430204e+01
6 7 8 9 10
4.401633e+01 -1.433905e+02 -9.819038e+00 -3.788190e+02 1.418952e+02
11 12 13 14 15
8.089525e+01 4.516353e-15 4.665897e+01 2.230875e+02 -1.599125e+02
16 17 18 19 20
3.329447e+02 2.212304e+02 -1.600553e+02 2.575379e+02 3.610932e+01
21 22 23 24 25
1.181093e+02 1.528236e+02 4.082361e+01 -3.090272e+02 1.215873e+02
26 27 28 29 30
7.015902e+00 -3.489841e+02 2.018730e+02 -1.284124e+01 -8.412696e+01
31 32 33 34 35
5.773138e+02 1.788852e+02 6.018852e+02 5.695995e+02 2.259951e+01
36 37 38 39 40
6.287487e+02 1.433632e+02 1.117918e+02 2.097918e+02 2.466489e+02
41 42 43 44 45
-4.110653e+02 3.864895e+01 -8.875785e+01 -4.221864e+02 -1.171864e+02
46 47 48 49 50
-5.284721e+02 -2.404721e+02 -1.873229e+02 -3.957084e+02 -1.742798e+02
51 52 53 54 55
-1.327984e+01 -2.964227e+02 -4.191370e+02 -2.574227e+02 -4.698295e+02
56 57 58 59 60
-2.712581e+02 -3.225806e+01 -5.535438e+02 -1.855438e+02 -3.939456e+01
61 62 63 64 65
-2.417800e+02 -2.383515e+02 1.686485e+02 -4.234943e+02 2.307914e+02
66 67 68 69 70
2.455057e+02 -2.589011e+02 1.546703e+02 1.266703e+02 4.038459e+01
71 72 73 74 75
3.163846e+02 8.453380e+01 4.641483e+02 1.557688e+01 3.905769e+02
76 77 78 79 80
4.243403e+01 3.767197e+02 1.734340e+02 1.260272e+02 3.335987e+02
81 82 83 84
-3.184013e+02 1.773129e+02 -3.468705e+01 -1.775378e+02
> postscript(file="/var/wessaorg/rcomp/tmp/6de5h1385408508.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 = 84
Frequency = 1
lag(myerror, k = 1) myerror
0 -1.382694e+02 NA
1 5.515918e+01 -1.382694e+02
2 -2.468408e+02 5.515918e+01
3 -1.039837e+02 -2.468408e+02
4 1.430204e+01 -1.039837e+02
5 4.401633e+01 1.430204e+01
6 -1.433905e+02 4.401633e+01
7 -9.819038e+00 -1.433905e+02
8 -3.788190e+02 -9.819038e+00
9 1.418952e+02 -3.788190e+02
10 8.089525e+01 1.418952e+02
11 4.516353e-15 8.089525e+01
12 4.665897e+01 4.516353e-15
13 2.230875e+02 4.665897e+01
14 -1.599125e+02 2.230875e+02
15 3.329447e+02 -1.599125e+02
16 2.212304e+02 3.329447e+02
17 -1.600553e+02 2.212304e+02
18 2.575379e+02 -1.600553e+02
19 3.610932e+01 2.575379e+02
20 1.181093e+02 3.610932e+01
21 1.528236e+02 1.181093e+02
22 4.082361e+01 1.528236e+02
23 -3.090272e+02 4.082361e+01
24 1.215873e+02 -3.090272e+02
25 7.015902e+00 1.215873e+02
26 -3.489841e+02 7.015902e+00
27 2.018730e+02 -3.489841e+02
28 -1.284124e+01 2.018730e+02
29 -8.412696e+01 -1.284124e+01
30 5.773138e+02 -8.412696e+01
31 1.788852e+02 5.773138e+02
32 6.018852e+02 1.788852e+02
33 5.695995e+02 6.018852e+02
34 2.259951e+01 5.695995e+02
35 6.287487e+02 2.259951e+01
36 1.433632e+02 6.287487e+02
37 1.117918e+02 1.433632e+02
38 2.097918e+02 1.117918e+02
39 2.466489e+02 2.097918e+02
40 -4.110653e+02 2.466489e+02
41 3.864895e+01 -4.110653e+02
42 -8.875785e+01 3.864895e+01
43 -4.221864e+02 -8.875785e+01
44 -1.171864e+02 -4.221864e+02
45 -5.284721e+02 -1.171864e+02
46 -2.404721e+02 -5.284721e+02
47 -1.873229e+02 -2.404721e+02
48 -3.957084e+02 -1.873229e+02
49 -1.742798e+02 -3.957084e+02
50 -1.327984e+01 -1.742798e+02
51 -2.964227e+02 -1.327984e+01
52 -4.191370e+02 -2.964227e+02
53 -2.574227e+02 -4.191370e+02
54 -4.698295e+02 -2.574227e+02
55 -2.712581e+02 -4.698295e+02
56 -3.225806e+01 -2.712581e+02
57 -5.535438e+02 -3.225806e+01
58 -1.855438e+02 -5.535438e+02
59 -3.939456e+01 -1.855438e+02
60 -2.417800e+02 -3.939456e+01
61 -2.383515e+02 -2.417800e+02
62 1.686485e+02 -2.383515e+02
63 -4.234943e+02 1.686485e+02
64 2.307914e+02 -4.234943e+02
65 2.455057e+02 2.307914e+02
66 -2.589011e+02 2.455057e+02
67 1.546703e+02 -2.589011e+02
68 1.266703e+02 1.546703e+02
69 4.038459e+01 1.266703e+02
70 3.163846e+02 4.038459e+01
71 8.453380e+01 3.163846e+02
72 4.641483e+02 8.453380e+01
73 1.557688e+01 4.641483e+02
74 3.905769e+02 1.557688e+01
75 4.243403e+01 3.905769e+02
76 3.767197e+02 4.243403e+01
77 1.734340e+02 3.767197e+02
78 1.260272e+02 1.734340e+02
79 3.335987e+02 1.260272e+02
80 -3.184013e+02 3.335987e+02
81 1.773129e+02 -3.184013e+02
82 -3.468705e+01 1.773129e+02
83 -1.775378e+02 -3.468705e+01
84 NA -1.775378e+02
> dum1 <- dum[2:length(myerror),]
> dum1
lag(myerror, k = 1) myerror
[1,] 5.515918e+01 -1.382694e+02
[2,] -2.468408e+02 5.515918e+01
[3,] -1.039837e+02 -2.468408e+02
[4,] 1.430204e+01 -1.039837e+02
[5,] 4.401633e+01 1.430204e+01
[6,] -1.433905e+02 4.401633e+01
[7,] -9.819038e+00 -1.433905e+02
[8,] -3.788190e+02 -9.819038e+00
[9,] 1.418952e+02 -3.788190e+02
[10,] 8.089525e+01 1.418952e+02
[11,] 4.516353e-15 8.089525e+01
[12,] 4.665897e+01 4.516353e-15
[13,] 2.230875e+02 4.665897e+01
[14,] -1.599125e+02 2.230875e+02
[15,] 3.329447e+02 -1.599125e+02
[16,] 2.212304e+02 3.329447e+02
[17,] -1.600553e+02 2.212304e+02
[18,] 2.575379e+02 -1.600553e+02
[19,] 3.610932e+01 2.575379e+02
[20,] 1.181093e+02 3.610932e+01
[21,] 1.528236e+02 1.181093e+02
[22,] 4.082361e+01 1.528236e+02
[23,] -3.090272e+02 4.082361e+01
[24,] 1.215873e+02 -3.090272e+02
[25,] 7.015902e+00 1.215873e+02
[26,] -3.489841e+02 7.015902e+00
[27,] 2.018730e+02 -3.489841e+02
[28,] -1.284124e+01 2.018730e+02
[29,] -8.412696e+01 -1.284124e+01
[30,] 5.773138e+02 -8.412696e+01
[31,] 1.788852e+02 5.773138e+02
[32,] 6.018852e+02 1.788852e+02
[33,] 5.695995e+02 6.018852e+02
[34,] 2.259951e+01 5.695995e+02
[35,] 6.287487e+02 2.259951e+01
[36,] 1.433632e+02 6.287487e+02
[37,] 1.117918e+02 1.433632e+02
[38,] 2.097918e+02 1.117918e+02
[39,] 2.466489e+02 2.097918e+02
[40,] -4.110653e+02 2.466489e+02
[41,] 3.864895e+01 -4.110653e+02
[42,] -8.875785e+01 3.864895e+01
[43,] -4.221864e+02 -8.875785e+01
[44,] -1.171864e+02 -4.221864e+02
[45,] -5.284721e+02 -1.171864e+02
[46,] -2.404721e+02 -5.284721e+02
[47,] -1.873229e+02 -2.404721e+02
[48,] -3.957084e+02 -1.873229e+02
[49,] -1.742798e+02 -3.957084e+02
[50,] -1.327984e+01 -1.742798e+02
[51,] -2.964227e+02 -1.327984e+01
[52,] -4.191370e+02 -2.964227e+02
[53,] -2.574227e+02 -4.191370e+02
[54,] -4.698295e+02 -2.574227e+02
[55,] -2.712581e+02 -4.698295e+02
[56,] -3.225806e+01 -2.712581e+02
[57,] -5.535438e+02 -3.225806e+01
[58,] -1.855438e+02 -5.535438e+02
[59,] -3.939456e+01 -1.855438e+02
[60,] -2.417800e+02 -3.939456e+01
[61,] -2.383515e+02 -2.417800e+02
[62,] 1.686485e+02 -2.383515e+02
[63,] -4.234943e+02 1.686485e+02
[64,] 2.307914e+02 -4.234943e+02
[65,] 2.455057e+02 2.307914e+02
[66,] -2.589011e+02 2.455057e+02
[67,] 1.546703e+02 -2.589011e+02
[68,] 1.266703e+02 1.546703e+02
[69,] 4.038459e+01 1.266703e+02
[70,] 3.163846e+02 4.038459e+01
[71,] 8.453380e+01 3.163846e+02
[72,] 4.641483e+02 8.453380e+01
[73,] 1.557688e+01 4.641483e+02
[74,] 3.905769e+02 1.557688e+01
[75,] 4.243403e+01 3.905769e+02
[76,] 3.767197e+02 4.243403e+01
[77,] 1.734340e+02 3.767197e+02
[78,] 1.260272e+02 1.734340e+02
[79,] 3.335987e+02 1.260272e+02
[80,] -3.184013e+02 3.335987e+02
[81,] 1.773129e+02 -3.184013e+02
[82,] -3.468705e+01 1.773129e+02
[83,] -1.775378e+02 -3.468705e+01
> z <- as.data.frame(dum1)
> z
lag(myerror, k = 1) myerror
1 5.515918e+01 -1.382694e+02
2 -2.468408e+02 5.515918e+01
3 -1.039837e+02 -2.468408e+02
4 1.430204e+01 -1.039837e+02
5 4.401633e+01 1.430204e+01
6 -1.433905e+02 4.401633e+01
7 -9.819038e+00 -1.433905e+02
8 -3.788190e+02 -9.819038e+00
9 1.418952e+02 -3.788190e+02
10 8.089525e+01 1.418952e+02
11 4.516353e-15 8.089525e+01
12 4.665897e+01 4.516353e-15
13 2.230875e+02 4.665897e+01
14 -1.599125e+02 2.230875e+02
15 3.329447e+02 -1.599125e+02
16 2.212304e+02 3.329447e+02
17 -1.600553e+02 2.212304e+02
18 2.575379e+02 -1.600553e+02
19 3.610932e+01 2.575379e+02
20 1.181093e+02 3.610932e+01
21 1.528236e+02 1.181093e+02
22 4.082361e+01 1.528236e+02
23 -3.090272e+02 4.082361e+01
24 1.215873e+02 -3.090272e+02
25 7.015902e+00 1.215873e+02
26 -3.489841e+02 7.015902e+00
27 2.018730e+02 -3.489841e+02
28 -1.284124e+01 2.018730e+02
29 -8.412696e+01 -1.284124e+01
30 5.773138e+02 -8.412696e+01
31 1.788852e+02 5.773138e+02
32 6.018852e+02 1.788852e+02
33 5.695995e+02 6.018852e+02
34 2.259951e+01 5.695995e+02
35 6.287487e+02 2.259951e+01
36 1.433632e+02 6.287487e+02
37 1.117918e+02 1.433632e+02
38 2.097918e+02 1.117918e+02
39 2.466489e+02 2.097918e+02
40 -4.110653e+02 2.466489e+02
41 3.864895e+01 -4.110653e+02
42 -8.875785e+01 3.864895e+01
43 -4.221864e+02 -8.875785e+01
44 -1.171864e+02 -4.221864e+02
45 -5.284721e+02 -1.171864e+02
46 -2.404721e+02 -5.284721e+02
47 -1.873229e+02 -2.404721e+02
48 -3.957084e+02 -1.873229e+02
49 -1.742798e+02 -3.957084e+02
50 -1.327984e+01 -1.742798e+02
51 -2.964227e+02 -1.327984e+01
52 -4.191370e+02 -2.964227e+02
53 -2.574227e+02 -4.191370e+02
54 -4.698295e+02 -2.574227e+02
55 -2.712581e+02 -4.698295e+02
56 -3.225806e+01 -2.712581e+02
57 -5.535438e+02 -3.225806e+01
58 -1.855438e+02 -5.535438e+02
59 -3.939456e+01 -1.855438e+02
60 -2.417800e+02 -3.939456e+01
61 -2.383515e+02 -2.417800e+02
62 1.686485e+02 -2.383515e+02
63 -4.234943e+02 1.686485e+02
64 2.307914e+02 -4.234943e+02
65 2.455057e+02 2.307914e+02
66 -2.589011e+02 2.455057e+02
67 1.546703e+02 -2.589011e+02
68 1.266703e+02 1.546703e+02
69 4.038459e+01 1.266703e+02
70 3.163846e+02 4.038459e+01
71 8.453380e+01 3.163846e+02
72 4.641483e+02 8.453380e+01
73 1.557688e+01 4.641483e+02
74 3.905769e+02 1.557688e+01
75 4.243403e+01 3.905769e+02
76 3.767197e+02 4.243403e+01
77 1.734340e+02 3.767197e+02
78 1.260272e+02 1.734340e+02
79 3.335987e+02 1.260272e+02
80 -3.184013e+02 3.335987e+02
81 1.773129e+02 -3.184013e+02
82 -3.468705e+01 1.773129e+02
83 -1.775378e+02 -3.468705e+01
> 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/7p9be1385408508.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/85o1e1385408508.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/9e7e21385408508.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')
Warning messages:
1: not plotting observations with leverage one:
12
2: not plotting observations with leverage one:
12
> par(opar)
> dev.off()
null device
1
> if (n > n25) {
+ postscript(file="/var/wessaorg/rcomp/tmp/1007um1385408508.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, 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/wessaorg/rcomp/tmp/11u28s1385408508.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/wessaorg/rcomp/tmp/124ro01385408508.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/wessaorg/rcomp/tmp/1343jk1385408508.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/wessaorg/rcomp/tmp/141p0i1385408509.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/wessaorg/rcomp/tmp/156ipx1385408509.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/wessaorg/rcomp/tmp/16uld01385408509.tab")
+ }
>
> try(system("convert tmp/1vhk41385408508.ps tmp/1vhk41385408508.png",intern=TRUE))
character(0)
> try(system("convert tmp/2pyn41385408508.ps tmp/2pyn41385408508.png",intern=TRUE))
character(0)
> try(system("convert tmp/3haag1385408508.ps tmp/3haag1385408508.png",intern=TRUE))
character(0)
> try(system("convert tmp/4td8r1385408508.ps tmp/4td8r1385408508.png",intern=TRUE))
character(0)
> try(system("convert tmp/56gy61385408508.ps tmp/56gy61385408508.png",intern=TRUE))
character(0)
> try(system("convert tmp/6de5h1385408508.ps tmp/6de5h1385408508.png",intern=TRUE))
character(0)
> try(system("convert tmp/7p9be1385408508.ps tmp/7p9be1385408508.png",intern=TRUE))
character(0)
> try(system("convert tmp/85o1e1385408508.ps tmp/85o1e1385408508.png",intern=TRUE))
character(0)
> try(system("convert tmp/9e7e21385408508.ps tmp/9e7e21385408508.png",intern=TRUE))
character(0)
> try(system("convert tmp/1007um1385408508.ps tmp/1007um1385408508.png",intern=TRUE))
character(0)
>
>
> proc.time()
user system elapsed
9.120 1.557 10.655