R version 2.12.0 (2010-10-15)
Copyright (C) 2010 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(4143,4429,5219,4929,5761,5592,4163,4962,5208,4755,4491,5732,5731,5040,6102,4904,5369,5578,4619,4731,5011,5299,4146,4625,4736,4219,5116,4205,4121,5103,4300,4578,3809,5657,4248,3830,4736,4839,4411,4570,4104,4801,3953,3828,4440,4026,4109,4785,3224,3552,3940,3913,3681,4309,3830,4143,4087,3818,3380,3430,3458,3970,5260,5024,5634,6549,4676),dim=c(1,67),dimnames=list(c('nb'),1:67))
> y <- array(NA,dim=c(1,67),dimnames=list(c('nb'),1:67))
> 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'
> #'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
nb M1 M2 M3 M4 M5 M6 M7 M8 M9 M10 M11 t
1 4143 1 0 0 0 0 0 0 0 0 0 0 1
2 4429 0 1 0 0 0 0 0 0 0 0 0 2
3 5219 0 0 1 0 0 0 0 0 0 0 0 3
4 4929 0 0 0 1 0 0 0 0 0 0 0 4
5 5761 0 0 0 0 1 0 0 0 0 0 0 5
6 5592 0 0 0 0 0 1 0 0 0 0 0 6
7 4163 0 0 0 0 0 0 1 0 0 0 0 7
8 4962 0 0 0 0 0 0 0 1 0 0 0 8
9 5208 0 0 0 0 0 0 0 0 1 0 0 9
10 4755 0 0 0 0 0 0 0 0 0 1 0 10
11 4491 0 0 0 0 0 0 0 0 0 0 1 11
12 5732 0 0 0 0 0 0 0 0 0 0 0 12
13 5731 1 0 0 0 0 0 0 0 0 0 0 13
14 5040 0 1 0 0 0 0 0 0 0 0 0 14
15 6102 0 0 1 0 0 0 0 0 0 0 0 15
16 4904 0 0 0 1 0 0 0 0 0 0 0 16
17 5369 0 0 0 0 1 0 0 0 0 0 0 17
18 5578 0 0 0 0 0 1 0 0 0 0 0 18
19 4619 0 0 0 0 0 0 1 0 0 0 0 19
20 4731 0 0 0 0 0 0 0 1 0 0 0 20
21 5011 0 0 0 0 0 0 0 0 1 0 0 21
22 5299 0 0 0 0 0 0 0 0 0 1 0 22
23 4146 0 0 0 0 0 0 0 0 0 0 1 23
24 4625 0 0 0 0 0 0 0 0 0 0 0 24
25 4736 1 0 0 0 0 0 0 0 0 0 0 25
26 4219 0 1 0 0 0 0 0 0 0 0 0 26
27 5116 0 0 1 0 0 0 0 0 0 0 0 27
28 4205 0 0 0 1 0 0 0 0 0 0 0 28
29 4121 0 0 0 0 1 0 0 0 0 0 0 29
30 5103 0 0 0 0 0 1 0 0 0 0 0 30
31 4300 0 0 0 0 0 0 1 0 0 0 0 31
32 4578 0 0 0 0 0 0 0 1 0 0 0 32
33 3809 0 0 0 0 0 0 0 0 1 0 0 33
34 5657 0 0 0 0 0 0 0 0 0 1 0 34
35 4248 0 0 0 0 0 0 0 0 0 0 1 35
36 3830 0 0 0 0 0 0 0 0 0 0 0 36
37 4736 1 0 0 0 0 0 0 0 0 0 0 37
38 4839 0 1 0 0 0 0 0 0 0 0 0 38
39 4411 0 0 1 0 0 0 0 0 0 0 0 39
40 4570 0 0 0 1 0 0 0 0 0 0 0 40
41 4104 0 0 0 0 1 0 0 0 0 0 0 41
42 4801 0 0 0 0 0 1 0 0 0 0 0 42
43 3953 0 0 0 0 0 0 1 0 0 0 0 43
44 3828 0 0 0 0 0 0 0 1 0 0 0 44
45 4440 0 0 0 0 0 0 0 0 1 0 0 45
46 4026 0 0 0 0 0 0 0 0 0 1 0 46
47 4109 0 0 0 0 0 0 0 0 0 0 1 47
48 4785 0 0 0 0 0 0 0 0 0 0 0 48
49 3224 1 0 0 0 0 0 0 0 0 0 0 49
50 3552 0 1 0 0 0 0 0 0 0 0 0 50
51 3940 0 0 1 0 0 0 0 0 0 0 0 51
52 3913 0 0 0 1 0 0 0 0 0 0 0 52
53 3681 0 0 0 0 1 0 0 0 0 0 0 53
54 4309 0 0 0 0 0 1 0 0 0 0 0 54
55 3830 0 0 0 0 0 0 1 0 0 0 0 55
56 4143 0 0 0 0 0 0 0 1 0 0 0 56
57 4087 0 0 0 0 0 0 0 0 1 0 0 57
58 3818 0 0 0 0 0 0 0 0 0 1 0 58
59 3380 0 0 0 0 0 0 0 0 0 0 1 59
60 3430 0 0 0 0 0 0 0 0 0 0 0 60
61 3458 1 0 0 0 0 0 0 0 0 0 0 61
62 3970 0 1 0 0 0 0 0 0 0 0 0 62
63 5260 0 0 1 0 0 0 0 0 0 0 0 63
64 5024 0 0 0 1 0 0 0 0 0 0 0 64
65 5634 0 0 0 0 1 0 0 0 0 0 0 65
66 6549 0 0 0 0 0 1 0 0 0 0 0 66
67 4676 0 0 0 0 0 0 1 0 0 0 0 67
> k <- length(x[1,])
> df <- as.data.frame(x)
> (mylm <- lm(df))
Call:
lm(formula = df)
Coefficients:
(Intercept) M1 M2 M3 M4 M5
5017.80 -217.04 -198.61 482.82 80.58 283.01
M6 M7 M8 M9 M10 M11
841.60 -208.64 -91.71 -14.18 200.74 -420.53
t
-14.93
> (mysum <- summary(mylm))
Call:
lm(formula = df)
Residuals:
Min 1Q Median 3Q Max
-845.30 -436.72 -34.53 321.42 1674.83
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 5017.800 309.380 16.219 < 2e-16 ***
M1 -217.039 373.571 -0.581 0.563668
M2 -198.611 373.387 -0.532 0.596965
M3 482.817 373.244 1.294 0.201320
M4 80.578 373.141 0.216 0.829845
M5 283.006 373.080 0.759 0.451411
M6 841.600 373.059 2.256 0.028147 *
M7 -208.639 373.080 -0.559 0.578314
M8 -91.711 389.961 -0.235 0.814959
M9 -14.183 389.824 -0.036 0.971110
M10 200.744 389.726 0.515 0.608592
M11 -420.528 389.667 -1.079 0.285297
t -14.928 3.909 -3.819 0.000348 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 616.1 on 54 degrees of freedom
Multiple R-squared: 0.3884, Adjusted R-squared: 0.2525
F-statistic: 2.858 on 12 and 54 DF, p-value: 0.004135
> 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.34921417 0.69842834 0.6507858
[2,] 0.42991873 0.85983746 0.5700813
[3,] 0.32714985 0.65429971 0.6728501
[4,] 0.20851398 0.41702797 0.7914860
[5,] 0.17063985 0.34127969 0.8293602
[6,] 0.13416069 0.26832138 0.8658393
[7,] 0.08803468 0.17606936 0.9119653
[8,] 0.07122954 0.14245908 0.9287705
[9,] 0.13852724 0.27705448 0.8614728
[10,] 0.11242964 0.22485928 0.8875704
[11,] 0.09270342 0.18540684 0.9072966
[12,] 0.07686269 0.15372537 0.9231373
[13,] 0.06399301 0.12798602 0.9360070
[14,] 0.10695249 0.21390497 0.8930475
[15,] 0.07123421 0.14246841 0.9287658
[16,] 0.04579260 0.09158520 0.9542074
[17,] 0.03005615 0.06011231 0.9699438
[18,] 0.03487936 0.06975873 0.9651206
[19,] 0.08175379 0.16350759 0.9182462
[20,] 0.06113805 0.12227610 0.9388620
[21,] 0.06504003 0.13008005 0.9349600
[22,] 0.10142921 0.20285841 0.8985708
[23,] 0.17062269 0.34124538 0.8293773
[24,] 0.14712085 0.29424169 0.8528792
[25,] 0.12356222 0.24712444 0.8764378
[26,] 0.09341023 0.18682045 0.9065898
[27,] 0.06050949 0.12101897 0.9394905
[28,] 0.03933898 0.07867797 0.9606610
[29,] 0.02541821 0.05083643 0.9745818
[30,] 0.02003786 0.04007573 0.9799621
[31,] 0.01875574 0.03751148 0.9812443
[32,] 0.03010512 0.06021024 0.9698949
[33,] 0.36227614 0.72455228 0.6377239
[34,] 0.49347990 0.98695980 0.5065201
[35,] 0.61995755 0.76008490 0.3800425
[36,] 0.46882239 0.93764478 0.5311776
> postscript(file="/var/www/rcomp/tmp/17wl51292882529.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/www/rcomp/tmp/2i5271292882529.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/www/rcomp/tmp/3i5271292882529.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/www/rcomp/tmp/4i5271292882529.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/www/rcomp/tmp/5i5271292882529.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 = 67
Frequency = 1
1 2 3 4 5 6 7
-642.83333 -360.33333 -236.83333 -109.66667 534.83333 -177.83333 -541.66667
8 9 10 11 12 13 14
155.33333 338.73333 -314.26667 57.93333 893.33333 1124.30000 429.80000
15 16 17 18 19 20 21
825.30000 44.46667 321.96667 -12.70000 93.46667 103.46667 320.86667
22 23 24 25 26 27 28
408.86667 -107.93333 -34.53333 308.43333 -212.06667 18.43333 -475.40000
29 30 31 32 33 34 35
-746.90000 -308.56667 -46.40000 129.60000 -702.00000 946.00000 173.20000
36 37 38 39 40 41 42
-650.40000 487.56667 587.06667 -507.43333 68.73333 -584.76667 -431.43333
43 44 45 46 47 48 49
-214.26667 -441.26667 108.13333 -505.86667 213.33333 483.73333 -845.30000
50 51 52 53 54 55 56
-520.80000 -799.30000 -409.13333 -828.63333 -744.30000 -158.13333 52.86667
57 58 59 60 61 62 63
-65.73333 -534.73333 -336.53333 -692.13333 -432.16667 76.33333 699.83333
64 65 66 67
881.00000 1303.50000 1674.83333 867.00000
> postscript(file="/var/www/rcomp/tmp/6xz9h1292882529.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 = 67
Frequency = 1
lag(myerror, k = 1) myerror
0 -642.83333 NA
1 -360.33333 -642.83333
2 -236.83333 -360.33333
3 -109.66667 -236.83333
4 534.83333 -109.66667
5 -177.83333 534.83333
6 -541.66667 -177.83333
7 155.33333 -541.66667
8 338.73333 155.33333
9 -314.26667 338.73333
10 57.93333 -314.26667
11 893.33333 57.93333
12 1124.30000 893.33333
13 429.80000 1124.30000
14 825.30000 429.80000
15 44.46667 825.30000
16 321.96667 44.46667
17 -12.70000 321.96667
18 93.46667 -12.70000
19 103.46667 93.46667
20 320.86667 103.46667
21 408.86667 320.86667
22 -107.93333 408.86667
23 -34.53333 -107.93333
24 308.43333 -34.53333
25 -212.06667 308.43333
26 18.43333 -212.06667
27 -475.40000 18.43333
28 -746.90000 -475.40000
29 -308.56667 -746.90000
30 -46.40000 -308.56667
31 129.60000 -46.40000
32 -702.00000 129.60000
33 946.00000 -702.00000
34 173.20000 946.00000
35 -650.40000 173.20000
36 487.56667 -650.40000
37 587.06667 487.56667
38 -507.43333 587.06667
39 68.73333 -507.43333
40 -584.76667 68.73333
41 -431.43333 -584.76667
42 -214.26667 -431.43333
43 -441.26667 -214.26667
44 108.13333 -441.26667
45 -505.86667 108.13333
46 213.33333 -505.86667
47 483.73333 213.33333
48 -845.30000 483.73333
49 -520.80000 -845.30000
50 -799.30000 -520.80000
51 -409.13333 -799.30000
52 -828.63333 -409.13333
53 -744.30000 -828.63333
54 -158.13333 -744.30000
55 52.86667 -158.13333
56 -65.73333 52.86667
57 -534.73333 -65.73333
58 -336.53333 -534.73333
59 -692.13333 -336.53333
60 -432.16667 -692.13333
61 76.33333 -432.16667
62 699.83333 76.33333
63 881.00000 699.83333
64 1303.50000 881.00000
65 1674.83333 1303.50000
66 867.00000 1674.83333
67 NA 867.00000
> dum1 <- dum[2:length(myerror),]
> dum1
lag(myerror, k = 1) myerror
[1,] -360.33333 -642.83333
[2,] -236.83333 -360.33333
[3,] -109.66667 -236.83333
[4,] 534.83333 -109.66667
[5,] -177.83333 534.83333
[6,] -541.66667 -177.83333
[7,] 155.33333 -541.66667
[8,] 338.73333 155.33333
[9,] -314.26667 338.73333
[10,] 57.93333 -314.26667
[11,] 893.33333 57.93333
[12,] 1124.30000 893.33333
[13,] 429.80000 1124.30000
[14,] 825.30000 429.80000
[15,] 44.46667 825.30000
[16,] 321.96667 44.46667
[17,] -12.70000 321.96667
[18,] 93.46667 -12.70000
[19,] 103.46667 93.46667
[20,] 320.86667 103.46667
[21,] 408.86667 320.86667
[22,] -107.93333 408.86667
[23,] -34.53333 -107.93333
[24,] 308.43333 -34.53333
[25,] -212.06667 308.43333
[26,] 18.43333 -212.06667
[27,] -475.40000 18.43333
[28,] -746.90000 -475.40000
[29,] -308.56667 -746.90000
[30,] -46.40000 -308.56667
[31,] 129.60000 -46.40000
[32,] -702.00000 129.60000
[33,] 946.00000 -702.00000
[34,] 173.20000 946.00000
[35,] -650.40000 173.20000
[36,] 487.56667 -650.40000
[37,] 587.06667 487.56667
[38,] -507.43333 587.06667
[39,] 68.73333 -507.43333
[40,] -584.76667 68.73333
[41,] -431.43333 -584.76667
[42,] -214.26667 -431.43333
[43,] -441.26667 -214.26667
[44,] 108.13333 -441.26667
[45,] -505.86667 108.13333
[46,] 213.33333 -505.86667
[47,] 483.73333 213.33333
[48,] -845.30000 483.73333
[49,] -520.80000 -845.30000
[50,] -799.30000 -520.80000
[51,] -409.13333 -799.30000
[52,] -828.63333 -409.13333
[53,] -744.30000 -828.63333
[54,] -158.13333 -744.30000
[55,] 52.86667 -158.13333
[56,] -65.73333 52.86667
[57,] -534.73333 -65.73333
[58,] -336.53333 -534.73333
[59,] -692.13333 -336.53333
[60,] -432.16667 -692.13333
[61,] 76.33333 -432.16667
[62,] 699.83333 76.33333
[63,] 881.00000 699.83333
[64,] 1303.50000 881.00000
[65,] 1674.83333 1303.50000
[66,] 867.00000 1674.83333
> z <- as.data.frame(dum1)
> z
lag(myerror, k = 1) myerror
1 -360.33333 -642.83333
2 -236.83333 -360.33333
3 -109.66667 -236.83333
4 534.83333 -109.66667
5 -177.83333 534.83333
6 -541.66667 -177.83333
7 155.33333 -541.66667
8 338.73333 155.33333
9 -314.26667 338.73333
10 57.93333 -314.26667
11 893.33333 57.93333
12 1124.30000 893.33333
13 429.80000 1124.30000
14 825.30000 429.80000
15 44.46667 825.30000
16 321.96667 44.46667
17 -12.70000 321.96667
18 93.46667 -12.70000
19 103.46667 93.46667
20 320.86667 103.46667
21 408.86667 320.86667
22 -107.93333 408.86667
23 -34.53333 -107.93333
24 308.43333 -34.53333
25 -212.06667 308.43333
26 18.43333 -212.06667
27 -475.40000 18.43333
28 -746.90000 -475.40000
29 -308.56667 -746.90000
30 -46.40000 -308.56667
31 129.60000 -46.40000
32 -702.00000 129.60000
33 946.00000 -702.00000
34 173.20000 946.00000
35 -650.40000 173.20000
36 487.56667 -650.40000
37 587.06667 487.56667
38 -507.43333 587.06667
39 68.73333 -507.43333
40 -584.76667 68.73333
41 -431.43333 -584.76667
42 -214.26667 -431.43333
43 -441.26667 -214.26667
44 108.13333 -441.26667
45 -505.86667 108.13333
46 213.33333 -505.86667
47 483.73333 213.33333
48 -845.30000 483.73333
49 -520.80000 -845.30000
50 -799.30000 -520.80000
51 -409.13333 -799.30000
52 -828.63333 -409.13333
53 -744.30000 -828.63333
54 -158.13333 -744.30000
55 52.86667 -158.13333
56 -65.73333 52.86667
57 -534.73333 -65.73333
58 -336.53333 -534.73333
59 -692.13333 -336.53333
60 -432.16667 -692.13333
61 76.33333 -432.16667
62 699.83333 76.33333
63 881.00000 699.83333
64 1303.50000 881.00000
65 1674.83333 1303.50000
66 867.00000 1674.83333
> plot(z,main=paste('Residual Lag plot, lowess, and regression line'), ylab='values of Residuals', xlab='lagged values of Residuals')
> lines(lowess(z))
> abline(lm(z))
> grid()
> dev.off()
null device
1
> postscript(file="/var/www/rcomp/tmp/7qrqk1292882529.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/www/rcomp/tmp/8qrqk1292882529.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/www/rcomp/tmp/910751292882529.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/www/rcomp/tmp/1010751292882529.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/www/rcomp/createtable file can be downloaded at http://www.wessa.net/cretab
> load(file="/var/www/rcomp/createtable")
>
> a<-table.start()
> a<-table.row.start(a)
> a<-table.element(a, 'Multiple Linear Regression - Estimated Regression Equation', 1, TRUE)
> a<-table.row.end(a)
> myeq <- colnames(x)[1]
> myeq <- paste(myeq, '[t] = ', sep='')
> for (i in 1:k){
+ if (mysum$coefficients[i,1] > 0) myeq <- paste(myeq, '+', '')
+ myeq <- paste(myeq, mysum$coefficients[i,1], sep=' ')
+ if (rownames(mysum$coefficients)[i] != '(Intercept)') {
+ myeq <- paste(myeq, rownames(mysum$coefficients)[i], sep='')
+ if (rownames(mysum$coefficients)[i] != 't') myeq <- paste(myeq, '[t]', sep='')
+ }
+ }
> myeq <- paste(myeq, ' + e[t]')
> a<-table.row.start(a)
> a<-table.element(a, myeq)
> a<-table.row.end(a)
> a<-table.end(a)
> table.save(a,file="/var/www/rcomp/tmp/11soyp1292882529.tab")
> a<-table.start()
> a<-table.row.start(a)
> a<-table.element(a,hyperlink('http://www.xycoon.com/ols1.htm','Multiple Linear Regression - Ordinary Least Squares',''), 6, TRUE)
> a<-table.row.end(a)
> a<-table.row.start(a)
> a<-table.element(a,'Variable',header=TRUE)
> a<-table.element(a,'Parameter',header=TRUE)
> a<-table.element(a,'S.D.',header=TRUE)
> a<-table.element(a,'T-STAT
H0: parameter = 0',header=TRUE)
> a<-table.element(a,'2-tail p-value',header=TRUE)
> a<-table.element(a,'1-tail p-value',header=TRUE)
> a<-table.row.end(a)
> for (i in 1:k){
+ a<-table.row.start(a)
+ a<-table.element(a,rownames(mysum$coefficients)[i],header=TRUE)
+ a<-table.element(a,mysum$coefficients[i,1])
+ a<-table.element(a, round(mysum$coefficients[i,2],6))
+ a<-table.element(a, round(mysum$coefficients[i,3],4))
+ a<-table.element(a, round(mysum$coefficients[i,4],6))
+ a<-table.element(a, round(mysum$coefficients[i,4]/2,6))
+ a<-table.row.end(a)
+ }
> a<-table.end(a)
> table.save(a,file="/var/www/rcomp/tmp/12lyfs1292882529.tab")
> a<-table.start()
> a<-table.row.start(a)
> a<-table.element(a, 'Multiple Linear Regression - Regression Statistics', 2, TRUE)
> a<-table.row.end(a)
> a<-table.row.start(a)
> a<-table.element(a, 'Multiple R',1,TRUE)
> a<-table.element(a, sqrt(mysum$r.squared))
> a<-table.row.end(a)
> a<-table.row.start(a)
> a<-table.element(a, 'R-squared',1,TRUE)
> a<-table.element(a, mysum$r.squared)
> a<-table.row.end(a)
> a<-table.row.start(a)
> a<-table.element(a, 'Adjusted R-squared',1,TRUE)
> a<-table.element(a, mysum$adj.r.squared)
> a<-table.row.end(a)
> a<-table.row.start(a)
> a<-table.element(a, 'F-TEST (value)',1,TRUE)
> a<-table.element(a, mysum$fstatistic[1])
> a<-table.row.end(a)
> a<-table.row.start(a)
> a<-table.element(a, 'F-TEST (DF numerator)',1,TRUE)
> a<-table.element(a, mysum$fstatistic[2])
> a<-table.row.end(a)
> a<-table.row.start(a)
> a<-table.element(a, 'F-TEST (DF denominator)',1,TRUE)
> a<-table.element(a, mysum$fstatistic[3])
> a<-table.row.end(a)
> a<-table.row.start(a)
> a<-table.element(a, 'p-value',1,TRUE)
> a<-table.element(a, 1-pf(mysum$fstatistic[1],mysum$fstatistic[2],mysum$fstatistic[3]))
> a<-table.row.end(a)
> a<-table.row.start(a)
> a<-table.element(a, 'Multiple Linear Regression - Residual Statistics', 2, TRUE)
> a<-table.row.end(a)
> a<-table.row.start(a)
> a<-table.element(a, 'Residual Standard Deviation',1,TRUE)
> a<-table.element(a, mysum$sigma)
> a<-table.row.end(a)
> a<-table.row.start(a)
> a<-table.element(a, 'Sum Squared Residuals',1,TRUE)
> a<-table.element(a, sum(myerror*myerror))
> a<-table.row.end(a)
> a<-table.end(a)
> table.save(a,file="/var/www/rcomp/tmp/139zum1292882529.tab")
> a<-table.start()
> a<-table.row.start(a)
> a<-table.element(a, 'Multiple Linear Regression - Actuals, Interpolation, and Residuals', 4, TRUE)
> a<-table.row.end(a)
> a<-table.row.start(a)
> a<-table.element(a, 'Time or Index', 1, TRUE)
> a<-table.element(a, 'Actuals', 1, TRUE)
> a<-table.element(a, 'Interpolation
Forecast', 1, TRUE)
> a<-table.element(a, 'Residuals
Prediction Error', 1, TRUE)
> a<-table.row.end(a)
> for (i in 1:n) {
+ a<-table.row.start(a)
+ a<-table.element(a,i, 1, TRUE)
+ a<-table.element(a,x[i])
+ a<-table.element(a,x[i]-mysum$resid[i])
+ a<-table.element(a,mysum$resid[i])
+ a<-table.row.end(a)
+ }
> a<-table.end(a)
> table.save(a,file="/var/www/rcomp/tmp/14k8t71292882529.tab")
> if (n > n25) {
+ a<-table.start()
+ a<-table.row.start(a)
+ a<-table.element(a,'Goldfeld-Quandt test for Heteroskedasticity',4,TRUE)
+ a<-table.row.end(a)
+ a<-table.row.start(a)
+ a<-table.element(a,'p-values',header=TRUE)
+ a<-table.element(a,'Alternative Hypothesis',3,header=TRUE)
+ a<-table.row.end(a)
+ a<-table.row.start(a)
+ a<-table.element(a,'breakpoint index',header=TRUE)
+ a<-table.element(a,'greater',header=TRUE)
+ a<-table.element(a,'2-sided',header=TRUE)
+ a<-table.element(a,'less',header=TRUE)
+ a<-table.row.end(a)
+ for (mypoint in kp3:nmkm3) {
+ a<-table.row.start(a)
+ a<-table.element(a,mypoint,header=TRUE)
+ a<-table.element(a,gqarr[mypoint-kp3+1,1])
+ a<-table.element(a,gqarr[mypoint-kp3+1,2])
+ a<-table.element(a,gqarr[mypoint-kp3+1,3])
+ a<-table.row.end(a)
+ }
+ a<-table.end(a)
+ table.save(a,file="/var/www/rcomp/tmp/15nrav1292882529.tab")
+ a<-table.start()
+ a<-table.row.start(a)
+ a<-table.element(a,'Meta Analysis of Goldfeld-Quandt test for Heteroskedasticity',4,TRUE)
+ a<-table.row.end(a)
+ a<-table.row.start(a)
+ a<-table.element(a,'Description',header=TRUE)
+ a<-table.element(a,'# significant tests',header=TRUE)
+ a<-table.element(a,'% significant tests',header=TRUE)
+ a<-table.element(a,'OK/NOK',header=TRUE)
+ a<-table.row.end(a)
+ a<-table.row.start(a)
+ a<-table.element(a,'1% type I error level',header=TRUE)
+ a<-table.element(a,numsignificant1)
+ a<-table.element(a,numsignificant1/numgqtests)
+ if (numsignificant1/numgqtests < 0.01) dum <- 'OK' else dum <- 'NOK'
+ a<-table.element(a,dum)
+ a<-table.row.end(a)
+ a<-table.row.start(a)
+ a<-table.element(a,'5% type I error level',header=TRUE)
+ a<-table.element(a,numsignificant5)
+ a<-table.element(a,numsignificant5/numgqtests)
+ if (numsignificant5/numgqtests < 0.05) dum <- 'OK' else dum <- 'NOK'
+ a<-table.element(a,dum)
+ a<-table.row.end(a)
+ a<-table.row.start(a)
+ a<-table.element(a,'10% type I error level',header=TRUE)
+ a<-table.element(a,numsignificant10)
+ a<-table.element(a,numsignificant10/numgqtests)
+ if (numsignificant10/numgqtests < 0.1) dum <- 'OK' else dum <- 'NOK'
+ a<-table.element(a,dum)
+ a<-table.row.end(a)
+ a<-table.end(a)
+ table.save(a,file="/var/www/rcomp/tmp/16k0841292882529.tab")
+ }
>
> try(system("convert tmp/17wl51292882529.ps tmp/17wl51292882529.png",intern=TRUE))
character(0)
> try(system("convert tmp/2i5271292882529.ps tmp/2i5271292882529.png",intern=TRUE))
character(0)
> try(system("convert tmp/3i5271292882529.ps tmp/3i5271292882529.png",intern=TRUE))
character(0)
> try(system("convert tmp/4i5271292882529.ps tmp/4i5271292882529.png",intern=TRUE))
character(0)
> try(system("convert tmp/5i5271292882529.ps tmp/5i5271292882529.png",intern=TRUE))
character(0)
> try(system("convert tmp/6xz9h1292882529.ps tmp/6xz9h1292882529.png",intern=TRUE))
character(0)
> try(system("convert tmp/7qrqk1292882529.ps tmp/7qrqk1292882529.png",intern=TRUE))
character(0)
> try(system("convert tmp/8qrqk1292882529.ps tmp/8qrqk1292882529.png",intern=TRUE))
character(0)
> try(system("convert tmp/910751292882529.ps tmp/910751292882529.png",intern=TRUE))
character(0)
> try(system("convert tmp/1010751292882529.ps tmp/1010751292882529.png",intern=TRUE))
character(0)
>
>
> proc.time()
user system elapsed
3.120 0.730 3.898