R version 2.9.0 (2009-04-17)
Copyright (C) 2009 The R Foundation for Statistical Computing
ISBN 3-900051-07-0
R is free software and comes with ABSOLUTELY NO WARRANTY.
You are welcome to redistribute it under certain conditions.
Type 'license()' or 'licence()' for distribution details.
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(8.4,420,8.4,418,8.4,410,8.6,418,8.9,426,8.8,428,8.3,430,7.5,424,7.2,423,7.4,427,8.8,441,9.3,449,9.3,452,8.7,462,8.2,455,8.3,461,8.5,461,8.6,463,8.5,462,8.2,456,8.1,455,7.9,456,8.6,472,8.7,472,8.7,471,8.5,465,8.4,459,8.5,465,8.7,468,8.7,467,8.6,463,8.5,460,8.3,462,8.00,461,8.2,476,8.1,476,8.1,471,8.00,453,7.9,443,7.9,442,8.00,444,8.00,438,7.9,427,8.00,424,7.7,416,7.2,406,7.5,431,7.3,434,7.00,418,7.00,412,7.00,404,7.2,409,7.3,412,7.1,406,6.8,398,6.4,397,6.1,385,6.5,390,7.7,413,7.9,413,7.5,401,6.9,397,6.6,397,6.9,409,7.7,419,8.00,424,8.00,428,7.7,430,7.3,424,7.4,433,8.1,456,8.3,459,8.2,446),dim=c(2,73),dimnames=list(c('wgb','nwwz'),1:73))
> y <- array(NA,dim=c(2,73),dimnames=list(c('wgb','nwwz'),1:73))
> 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
Attaching package: 'zoo'
The following object(s) are masked from package:base :
as.Date.numeric
> n25 <- 25 #minimum number of obs. for Goldfeld-Quandt test
> par1 <- as.numeric(par1)
> x <- t(y)
> k <- length(x[1,])
> n <- length(x[,1])
> x1 <- cbind(x[,par1], x[,1:k!=par1])
> mycolnames <- c(colnames(x)[par1], colnames(x)[1:k!=par1])
> colnames(x1) <- mycolnames #colnames(x)[par1]
> x <- x1
> if (par3 == 'First Differences'){
+ x2 <- array(0, dim=c(n-1,k), dimnames=list(1:(n-1), paste('(1-B)',colnames(x),sep='')))
+ for (i in 1:n-1) {
+ for (j in 1:k) {
+ x2[i,j] <- x[i+1,j] - x[i,j]
+ }
+ }
+ x <- x2
+ }
> if (par2 == 'Include Monthly Dummies'){
+ x2 <- array(0, dim=c(n,11), dimnames=list(1:n, paste('M', seq(1:11), sep ='')))
+ for (i in 1:11){
+ x2[seq(i,n,12),i] <- 1
+ }
+ x <- cbind(x, x2)
+ }
> if (par2 == 'Include Quarterly Dummies'){
+ x2 <- array(0, dim=c(n,3), dimnames=list(1:n, paste('Q', seq(1:3), sep ='')))
+ for (i in 1:3){
+ x2[seq(i,n,4),i] <- 1
+ }
+ x <- cbind(x, x2)
+ }
> k <- length(x[1,])
> if (par3 == 'Linear Trend'){
+ x <- cbind(x, c(1:n))
+ colnames(x)[k+1] <- 't'
+ }
> x
wgb nwwz M1 M2 M3 M4 M5 M6 M7 M8 M9 M10 M11 t
1 8.4 420 1 0 0 0 0 0 0 0 0 0 0 1
2 8.4 418 0 1 0 0 0 0 0 0 0 0 0 2
3 8.4 410 0 0 1 0 0 0 0 0 0 0 0 3
4 8.6 418 0 0 0 1 0 0 0 0 0 0 0 4
5 8.9 426 0 0 0 0 1 0 0 0 0 0 0 5
6 8.8 428 0 0 0 0 0 1 0 0 0 0 0 6
7 8.3 430 0 0 0 0 0 0 1 0 0 0 0 7
8 7.5 424 0 0 0 0 0 0 0 1 0 0 0 8
9 7.2 423 0 0 0 0 0 0 0 0 1 0 0 9
10 7.4 427 0 0 0 0 0 0 0 0 0 1 0 10
11 8.8 441 0 0 0 0 0 0 0 0 0 0 1 11
12 9.3 449 0 0 0 0 0 0 0 0 0 0 0 12
13 9.3 452 1 0 0 0 0 0 0 0 0 0 0 13
14 8.7 462 0 1 0 0 0 0 0 0 0 0 0 14
15 8.2 455 0 0 1 0 0 0 0 0 0 0 0 15
16 8.3 461 0 0 0 1 0 0 0 0 0 0 0 16
17 8.5 461 0 0 0 0 1 0 0 0 0 0 0 17
18 8.6 463 0 0 0 0 0 1 0 0 0 0 0 18
19 8.5 462 0 0 0 0 0 0 1 0 0 0 0 19
20 8.2 456 0 0 0 0 0 0 0 1 0 0 0 20
21 8.1 455 0 0 0 0 0 0 0 0 1 0 0 21
22 7.9 456 0 0 0 0 0 0 0 0 0 1 0 22
23 8.6 472 0 0 0 0 0 0 0 0 0 0 1 23
24 8.7 472 0 0 0 0 0 0 0 0 0 0 0 24
25 8.7 471 1 0 0 0 0 0 0 0 0 0 0 25
26 8.5 465 0 1 0 0 0 0 0 0 0 0 0 26
27 8.4 459 0 0 1 0 0 0 0 0 0 0 0 27
28 8.5 465 0 0 0 1 0 0 0 0 0 0 0 28
29 8.7 468 0 0 0 0 1 0 0 0 0 0 0 29
30 8.7 467 0 0 0 0 0 1 0 0 0 0 0 30
31 8.6 463 0 0 0 0 0 0 1 0 0 0 0 31
32 8.5 460 0 0 0 0 0 0 0 1 0 0 0 32
33 8.3 462 0 0 0 0 0 0 0 0 1 0 0 33
34 8.0 461 0 0 0 0 0 0 0 0 0 1 0 34
35 8.2 476 0 0 0 0 0 0 0 0 0 0 1 35
36 8.1 476 0 0 0 0 0 0 0 0 0 0 0 36
37 8.1 471 1 0 0 0 0 0 0 0 0 0 0 37
38 8.0 453 0 1 0 0 0 0 0 0 0 0 0 38
39 7.9 443 0 0 1 0 0 0 0 0 0 0 0 39
40 7.9 442 0 0 0 1 0 0 0 0 0 0 0 40
41 8.0 444 0 0 0 0 1 0 0 0 0 0 0 41
42 8.0 438 0 0 0 0 0 1 0 0 0 0 0 42
43 7.9 427 0 0 0 0 0 0 1 0 0 0 0 43
44 8.0 424 0 0 0 0 0 0 0 1 0 0 0 44
45 7.7 416 0 0 0 0 0 0 0 0 1 0 0 45
46 7.2 406 0 0 0 0 0 0 0 0 0 1 0 46
47 7.5 431 0 0 0 0 0 0 0 0 0 0 1 47
48 7.3 434 0 0 0 0 0 0 0 0 0 0 0 48
49 7.0 418 1 0 0 0 0 0 0 0 0 0 0 49
50 7.0 412 0 1 0 0 0 0 0 0 0 0 0 50
51 7.0 404 0 0 1 0 0 0 0 0 0 0 0 51
52 7.2 409 0 0 0 1 0 0 0 0 0 0 0 52
53 7.3 412 0 0 0 0 1 0 0 0 0 0 0 53
54 7.1 406 0 0 0 0 0 1 0 0 0 0 0 54
55 6.8 398 0 0 0 0 0 0 1 0 0 0 0 55
56 6.4 397 0 0 0 0 0 0 0 1 0 0 0 56
57 6.1 385 0 0 0 0 0 0 0 0 1 0 0 57
58 6.5 390 0 0 0 0 0 0 0 0 0 1 0 58
59 7.7 413 0 0 0 0 0 0 0 0 0 0 1 59
60 7.9 413 0 0 0 0 0 0 0 0 0 0 0 60
61 7.5 401 1 0 0 0 0 0 0 0 0 0 0 61
62 6.9 397 0 1 0 0 0 0 0 0 0 0 0 62
63 6.6 397 0 0 1 0 0 0 0 0 0 0 0 63
64 6.9 409 0 0 0 1 0 0 0 0 0 0 0 64
65 7.7 419 0 0 0 0 1 0 0 0 0 0 0 65
66 8.0 424 0 0 0 0 0 1 0 0 0 0 0 66
67 8.0 428 0 0 0 0 0 0 1 0 0 0 0 67
68 7.7 430 0 0 0 0 0 0 0 1 0 0 0 68
69 7.3 424 0 0 0 0 0 0 0 0 1 0 0 69
70 7.4 433 0 0 0 0 0 0 0 0 0 1 0 70
71 8.1 456 0 0 0 0 0 0 0 0 0 0 1 71
72 8.3 459 0 0 0 0 0 0 0 0 0 0 0 72
73 8.2 446 1 0 0 0 0 0 0 0 0 0 0 73
> k <- length(x[1,])
> df <- as.data.frame(x)
> (mylm <- lm(df))
Call:
lm(formula = df)
Coefficients:
(Intercept) nwwz M1 M2 M3 M4
2.163900 0.014799 -0.004884 -0.247512 -0.304559 -0.229922
M5 M6 M7 M8 M9 M10
0.002713 0.042675 -0.082834 -0.327475 -0.516585 -0.572888
M11 t
-0.095566 -0.013429
> (mysum <- summary(mylm))
Call:
lm(formula = df)
Residuals:
Min 1Q Median 3Q Max
-0.68682 -0.22716 -0.01467 0.30467 0.65266
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 2.163900 0.987587 2.191 0.0324 *
nwwz 0.014799 0.002065 7.167 1.41e-09 ***
M1 -0.004884 0.211738 -0.023 0.9817
M2 -0.247512 0.222861 -1.111 0.2712
M3 -0.304559 0.225395 -1.351 0.1818
M4 -0.229922 0.222321 -1.034 0.3053
M5 0.002713 0.220512 0.012 0.9902
M6 0.042675 0.220438 0.194 0.8472
M7 -0.082834 0.221153 -0.375 0.7093
M8 -0.327475 0.221942 -1.476 0.1454
M9 -0.516585 0.223533 -2.311 0.0243 *
M10 -0.572888 0.222669 -2.573 0.0126 *
M11 -0.095566 0.217763 -0.439 0.6624
t -0.013429 0.002333 -5.757 3.26e-07 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 0.377 on 59 degrees of freedom
Multiple R-squared: 0.7588, Adjusted R-squared: 0.7057
F-statistic: 14.28 on 13 and 59 DF, p-value: 1.040e-13
> 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.8325751 0.3348497 0.16742485
[2,] 0.7843826 0.4312347 0.21561737
[3,] 0.6656708 0.6686584 0.33432920
[4,] 0.6298818 0.7402365 0.37011823
[5,] 0.6304104 0.7391791 0.36958957
[6,] 0.5211754 0.9576493 0.47882463
[7,] 0.5783578 0.8432845 0.42164223
[8,] 0.7997502 0.4004996 0.20024982
[9,] 0.7644482 0.4711037 0.23555183
[10,] 0.6963362 0.6073277 0.30366384
[11,] 0.6281449 0.7437102 0.37185510
[12,] 0.5492020 0.9015960 0.45079802
[13,] 0.4638890 0.9277781 0.53611095
[14,] 0.3782202 0.7564404 0.62177980
[15,] 0.3111837 0.6223674 0.68881628
[16,] 0.3486052 0.6972105 0.65139477
[17,] 0.3513113 0.7026226 0.64868868
[18,] 0.2794573 0.5589146 0.72054268
[19,] 0.3442425 0.6884849 0.65575754
[20,] 0.5472239 0.9055522 0.45277611
[21,] 0.6341639 0.7316722 0.36583610
[22,] 0.5638098 0.8723805 0.43619025
[23,] 0.4806830 0.9613659 0.51931705
[24,] 0.3985767 0.7971534 0.60142330
[25,] 0.3283194 0.6566389 0.67168056
[26,] 0.2565017 0.5130033 0.74349834
[27,] 0.1992483 0.3984966 0.80075169
[28,] 0.2982408 0.5964816 0.70175921
[29,] 0.5454162 0.9091676 0.45458381
[30,] 0.6565316 0.6869368 0.34346841
[31,] 0.5993345 0.8013310 0.40066550
[32,] 0.6285734 0.7428531 0.37142656
[33,] 0.6949246 0.6101509 0.30507544
[34,] 0.6126579 0.7746842 0.38734209
[35,] 0.5824999 0.8350002 0.41750009
[36,] 0.7466176 0.5067648 0.25338241
[37,] 0.8510292 0.2979417 0.14897084
[38,] 0.9460391 0.1079219 0.05396093
[39,] 0.9239787 0.1520426 0.07602131
[40,] 0.8696583 0.2606833 0.13034166
> postscript(file="/var/www/html/rcomp/tmp/1pgqd1258624337.ps",horizontal=F,pagecentre=F,paper="special",width=8.3333333333333,height=5.5555555555556)
> plot(x[,1], type='l', main='Actuals and Interpolation', ylab='value of Actuals and Interpolation (dots)', xlab='time or index')
> points(x[,1]-mysum$resid)
> grid()
> dev.off()
null device
1
> postscript(file="/var/www/html/rcomp/tmp/2h00c1258624337.ps",horizontal=F,pagecentre=F,paper="special",width=8.3333333333333,height=5.5555555555556)
> plot(mysum$resid, type='b', pch=19, main='Residuals', ylab='value of Residuals', xlab='time or index')
> grid()
> dev.off()
null device
1
> postscript(file="/var/www/html/rcomp/tmp/3fpcv1258624337.ps",horizontal=F,pagecentre=F,paper="special",width=8.3333333333333,height=5.5555555555556)
> hist(mysum$resid, main='Residual Histogram', xlab='values of Residuals')
> grid()
> dev.off()
null device
1
> postscript(file="/var/www/html/rcomp/tmp/4sodh1258624337.ps",horizontal=F,pagecentre=F,paper="special",width=8.3333333333333,height=5.5555555555556)
> densityplot(~mysum$resid,col='black',main='Residual Density Plot', xlab='values of Residuals')
> dev.off()
null device
1
> postscript(file="/var/www/html/rcomp/tmp/5f7zb1258624337.ps",horizontal=F,pagecentre=F,paper="special",width=8.3333333333333,height=5.5555555555556)
> qqnorm(mysum$resid, main='Residual Normal Q-Q Plot')
> qqline(mysum$resid)
> grid()
> dev.off()
null device
1
> (myerror <- as.ts(mysum$resid))
Time Series:
Start = 1
End = 73
Frequency = 1
1 2 3 4 5 6
0.038983714 0.324638874 0.513503500 0.533906221 0.496311209 0.340181504
7 8 9 10 11 12
-0.050478361 -0.503616002 -0.586278134 -0.375741173 0.353184905 0.652659281
13 14 15 16 17 18
0.626576053 0.134647539 -0.191286475 -0.241286475 -0.260492370 -0.216622076
19 20 21 22 23 24
-0.162886021 -0.116023663 0.001314205 -0.143752915 -0.144424116 -0.126560624
25 26 27 28 29 30
-0.093449293 0.051400425 0.110667772 0.060667772 -0.002934042 -0.014667829
31 32 33 34 35 36
0.083464144 0.285930584 0.258872534 0.043402692 -0.442469870 -0.624606377
37 38 39 40 41 42
-0.532300488 -0.109867096 0.008594809 -0.037814714 -0.186617888 -0.124358477
43 44 45 46 47 48
0.077363972 0.479830412 0.500758757 0.218476672 -0.315382286 -0.641914712
49 50 51 52 53 54
-0.686823788 -0.341974070 -0.153109444 -0.088310804 -0.251912618 -0.389653207
55 56 57 58 59 60
-0.432326676 -0.559457516 -0.479334612 -0.083596291 0.312142031 0.430005523
61 62 63 64 65 66
0.225901889 -0.058845672 -0.288370162 -0.227162000 0.205645710 0.405120086
67 68 69 70 71 72
0.484862942 0.413336185 0.304667250 0.341211014 0.236949336 0.310416909
73
0.421111914
> postscript(file="/var/www/html/rcomp/tmp/64j3g1258624337.ps",horizontal=F,pagecentre=F,paper="special",width=8.3333333333333,height=5.5555555555556)
> dum <- cbind(lag(myerror,k=1),myerror)
> dum
Time Series:
Start = 0
End = 73
Frequency = 1
lag(myerror, k = 1) myerror
0 0.038983714 NA
1 0.324638874 0.038983714
2 0.513503500 0.324638874
3 0.533906221 0.513503500
4 0.496311209 0.533906221
5 0.340181504 0.496311209
6 -0.050478361 0.340181504
7 -0.503616002 -0.050478361
8 -0.586278134 -0.503616002
9 -0.375741173 -0.586278134
10 0.353184905 -0.375741173
11 0.652659281 0.353184905
12 0.626576053 0.652659281
13 0.134647539 0.626576053
14 -0.191286475 0.134647539
15 -0.241286475 -0.191286475
16 -0.260492370 -0.241286475
17 -0.216622076 -0.260492370
18 -0.162886021 -0.216622076
19 -0.116023663 -0.162886021
20 0.001314205 -0.116023663
21 -0.143752915 0.001314205
22 -0.144424116 -0.143752915
23 -0.126560624 -0.144424116
24 -0.093449293 -0.126560624
25 0.051400425 -0.093449293
26 0.110667772 0.051400425
27 0.060667772 0.110667772
28 -0.002934042 0.060667772
29 -0.014667829 -0.002934042
30 0.083464144 -0.014667829
31 0.285930584 0.083464144
32 0.258872534 0.285930584
33 0.043402692 0.258872534
34 -0.442469870 0.043402692
35 -0.624606377 -0.442469870
36 -0.532300488 -0.624606377
37 -0.109867096 -0.532300488
38 0.008594809 -0.109867096
39 -0.037814714 0.008594809
40 -0.186617888 -0.037814714
41 -0.124358477 -0.186617888
42 0.077363972 -0.124358477
43 0.479830412 0.077363972
44 0.500758757 0.479830412
45 0.218476672 0.500758757
46 -0.315382286 0.218476672
47 -0.641914712 -0.315382286
48 -0.686823788 -0.641914712
49 -0.341974070 -0.686823788
50 -0.153109444 -0.341974070
51 -0.088310804 -0.153109444
52 -0.251912618 -0.088310804
53 -0.389653207 -0.251912618
54 -0.432326676 -0.389653207
55 -0.559457516 -0.432326676
56 -0.479334612 -0.559457516
57 -0.083596291 -0.479334612
58 0.312142031 -0.083596291
59 0.430005523 0.312142031
60 0.225901889 0.430005523
61 -0.058845672 0.225901889
62 -0.288370162 -0.058845672
63 -0.227162000 -0.288370162
64 0.205645710 -0.227162000
65 0.405120086 0.205645710
66 0.484862942 0.405120086
67 0.413336185 0.484862942
68 0.304667250 0.413336185
69 0.341211014 0.304667250
70 0.236949336 0.341211014
71 0.310416909 0.236949336
72 0.421111914 0.310416909
73 NA 0.421111914
> dum1 <- dum[2:length(myerror),]
> dum1
lag(myerror, k = 1) myerror
[1,] 0.324638874 0.038983714
[2,] 0.513503500 0.324638874
[3,] 0.533906221 0.513503500
[4,] 0.496311209 0.533906221
[5,] 0.340181504 0.496311209
[6,] -0.050478361 0.340181504
[7,] -0.503616002 -0.050478361
[8,] -0.586278134 -0.503616002
[9,] -0.375741173 -0.586278134
[10,] 0.353184905 -0.375741173
[11,] 0.652659281 0.353184905
[12,] 0.626576053 0.652659281
[13,] 0.134647539 0.626576053
[14,] -0.191286475 0.134647539
[15,] -0.241286475 -0.191286475
[16,] -0.260492370 -0.241286475
[17,] -0.216622076 -0.260492370
[18,] -0.162886021 -0.216622076
[19,] -0.116023663 -0.162886021
[20,] 0.001314205 -0.116023663
[21,] -0.143752915 0.001314205
[22,] -0.144424116 -0.143752915
[23,] -0.126560624 -0.144424116
[24,] -0.093449293 -0.126560624
[25,] 0.051400425 -0.093449293
[26,] 0.110667772 0.051400425
[27,] 0.060667772 0.110667772
[28,] -0.002934042 0.060667772
[29,] -0.014667829 -0.002934042
[30,] 0.083464144 -0.014667829
[31,] 0.285930584 0.083464144
[32,] 0.258872534 0.285930584
[33,] 0.043402692 0.258872534
[34,] -0.442469870 0.043402692
[35,] -0.624606377 -0.442469870
[36,] -0.532300488 -0.624606377
[37,] -0.109867096 -0.532300488
[38,] 0.008594809 -0.109867096
[39,] -0.037814714 0.008594809
[40,] -0.186617888 -0.037814714
[41,] -0.124358477 -0.186617888
[42,] 0.077363972 -0.124358477
[43,] 0.479830412 0.077363972
[44,] 0.500758757 0.479830412
[45,] 0.218476672 0.500758757
[46,] -0.315382286 0.218476672
[47,] -0.641914712 -0.315382286
[48,] -0.686823788 -0.641914712
[49,] -0.341974070 -0.686823788
[50,] -0.153109444 -0.341974070
[51,] -0.088310804 -0.153109444
[52,] -0.251912618 -0.088310804
[53,] -0.389653207 -0.251912618
[54,] -0.432326676 -0.389653207
[55,] -0.559457516 -0.432326676
[56,] -0.479334612 -0.559457516
[57,] -0.083596291 -0.479334612
[58,] 0.312142031 -0.083596291
[59,] 0.430005523 0.312142031
[60,] 0.225901889 0.430005523
[61,] -0.058845672 0.225901889
[62,] -0.288370162 -0.058845672
[63,] -0.227162000 -0.288370162
[64,] 0.205645710 -0.227162000
[65,] 0.405120086 0.205645710
[66,] 0.484862942 0.405120086
[67,] 0.413336185 0.484862942
[68,] 0.304667250 0.413336185
[69,] 0.341211014 0.304667250
[70,] 0.236949336 0.341211014
[71,] 0.310416909 0.236949336
[72,] 0.421111914 0.310416909
> z <- as.data.frame(dum1)
> z
lag(myerror, k = 1) myerror
1 0.324638874 0.038983714
2 0.513503500 0.324638874
3 0.533906221 0.513503500
4 0.496311209 0.533906221
5 0.340181504 0.496311209
6 -0.050478361 0.340181504
7 -0.503616002 -0.050478361
8 -0.586278134 -0.503616002
9 -0.375741173 -0.586278134
10 0.353184905 -0.375741173
11 0.652659281 0.353184905
12 0.626576053 0.652659281
13 0.134647539 0.626576053
14 -0.191286475 0.134647539
15 -0.241286475 -0.191286475
16 -0.260492370 -0.241286475
17 -0.216622076 -0.260492370
18 -0.162886021 -0.216622076
19 -0.116023663 -0.162886021
20 0.001314205 -0.116023663
21 -0.143752915 0.001314205
22 -0.144424116 -0.143752915
23 -0.126560624 -0.144424116
24 -0.093449293 -0.126560624
25 0.051400425 -0.093449293
26 0.110667772 0.051400425
27 0.060667772 0.110667772
28 -0.002934042 0.060667772
29 -0.014667829 -0.002934042
30 0.083464144 -0.014667829
31 0.285930584 0.083464144
32 0.258872534 0.285930584
33 0.043402692 0.258872534
34 -0.442469870 0.043402692
35 -0.624606377 -0.442469870
36 -0.532300488 -0.624606377
37 -0.109867096 -0.532300488
38 0.008594809 -0.109867096
39 -0.037814714 0.008594809
40 -0.186617888 -0.037814714
41 -0.124358477 -0.186617888
42 0.077363972 -0.124358477
43 0.479830412 0.077363972
44 0.500758757 0.479830412
45 0.218476672 0.500758757
46 -0.315382286 0.218476672
47 -0.641914712 -0.315382286
48 -0.686823788 -0.641914712
49 -0.341974070 -0.686823788
50 -0.153109444 -0.341974070
51 -0.088310804 -0.153109444
52 -0.251912618 -0.088310804
53 -0.389653207 -0.251912618
54 -0.432326676 -0.389653207
55 -0.559457516 -0.432326676
56 -0.479334612 -0.559457516
57 -0.083596291 -0.479334612
58 0.312142031 -0.083596291
59 0.430005523 0.312142031
60 0.225901889 0.430005523
61 -0.058845672 0.225901889
62 -0.288370162 -0.058845672
63 -0.227162000 -0.288370162
64 0.205645710 -0.227162000
65 0.405120086 0.205645710
66 0.484862942 0.405120086
67 0.413336185 0.484862942
68 0.304667250 0.413336185
69 0.341211014 0.304667250
70 0.236949336 0.341211014
71 0.310416909 0.236949336
72 0.421111914 0.310416909
> plot(z,main=paste('Residual Lag plot, lowess, and regression line'), ylab='values of Residuals', xlab='lagged values of Residuals')
> lines(lowess(z))
> abline(lm(z))
> grid()
> dev.off()
null device
1
> postscript(file="/var/www/html/rcomp/tmp/7lsm51258624337.ps",horizontal=F,pagecentre=F,paper="special",width=8.3333333333333,height=5.5555555555556)
> acf(mysum$resid, lag.max=length(mysum$resid)/2, main='Residual Autocorrelation Function')
> grid()
> dev.off()
null device
1
> postscript(file="/var/www/html/rcomp/tmp/8uh5v1258624337.ps",horizontal=F,pagecentre=F,paper="special",width=8.3333333333333,height=5.5555555555556)
> pacf(mysum$resid, lag.max=length(mysum$resid)/2, main='Residual Partial Autocorrelation Function')
> grid()
> dev.off()
null device
1
> postscript(file="/var/www/html/rcomp/tmp/9395n1258624337.ps",horizontal=F,pagecentre=F,paper="special",width=8.3333333333333,height=5.5555555555556)
> opar <- par(mfrow = c(2,2), oma = c(0, 0, 1.1, 0))
> plot(mylm, las = 1, sub='Residual Diagnostics')
> par(opar)
> dev.off()
null device
1
> if (n > n25) {
+ postscript(file="/var/www/html/rcomp/tmp/10ealw1258624337.ps",horizontal=F,pagecentre=F,paper="special",width=8.3333333333333,height=5.5555555555556)
+ plot(kp3:nmkm3,gqarr[,2], main='Goldfeld-Quandt test',ylab='2-sided p-value',xlab='breakpoint')
+ grid()
+ dev.off()
+ }
null device
1
>
> #Note: the /var/www/html/rcomp/createtable file can be downloaded at http://www.wessa.net/cretab
> load(file="/var/www/html/rcomp/createtable")
>
> a<-table.start()
> a<-table.row.start(a)
> a<-table.element(a, 'Multiple Linear Regression - Estimated Regression Equation', 1, TRUE)
> a<-table.row.end(a)
> myeq <- colnames(x)[1]
> myeq <- paste(myeq, '[t] = ', sep='')
> for (i in 1:k){
+ if (mysum$coefficients[i,1] > 0) myeq <- paste(myeq, '+', '')
+ myeq <- paste(myeq, mysum$coefficients[i,1], sep=' ')
+ if (rownames(mysum$coefficients)[i] != '(Intercept)') {
+ myeq <- paste(myeq, rownames(mysum$coefficients)[i], sep='')
+ if (rownames(mysum$coefficients)[i] != 't') myeq <- paste(myeq, '[t]', sep='')
+ }
+ }
> myeq <- paste(myeq, ' + e[t]')
> a<-table.row.start(a)
> a<-table.element(a, myeq)
> a<-table.row.end(a)
> a<-table.end(a)
> table.save(a,file="/var/www/html/rcomp/tmp/11yx131258624337.tab")
> a<-table.start()
> a<-table.row.start(a)
> a<-table.element(a,hyperlink('http://www.xycoon.com/ols1.htm','Multiple Linear Regression - Ordinary Least Squares',''), 6, TRUE)
> a<-table.row.end(a)
> a<-table.row.start(a)
> a<-table.element(a,'Variable',header=TRUE)
> a<-table.element(a,'Parameter',header=TRUE)
> a<-table.element(a,'S.D.',header=TRUE)
> a<-table.element(a,'T-STAT
H0: parameter = 0',header=TRUE)
> a<-table.element(a,'2-tail p-value',header=TRUE)
> a<-table.element(a,'1-tail p-value',header=TRUE)
> a<-table.row.end(a)
> for (i in 1:k){
+ a<-table.row.start(a)
+ a<-table.element(a,rownames(mysum$coefficients)[i],header=TRUE)
+ a<-table.element(a,mysum$coefficients[i,1])
+ a<-table.element(a, round(mysum$coefficients[i,2],6))
+ a<-table.element(a, round(mysum$coefficients[i,3],4))
+ a<-table.element(a, round(mysum$coefficients[i,4],6))
+ a<-table.element(a, round(mysum$coefficients[i,4]/2,6))
+ a<-table.row.end(a)
+ }
> a<-table.end(a)
> table.save(a,file="/var/www/html/rcomp/tmp/12ya7m1258624337.tab")
> a<-table.start()
> a<-table.row.start(a)
> a<-table.element(a, 'Multiple Linear Regression - Regression Statistics', 2, TRUE)
> a<-table.row.end(a)
> a<-table.row.start(a)
> a<-table.element(a, 'Multiple R',1,TRUE)
> a<-table.element(a, sqrt(mysum$r.squared))
> a<-table.row.end(a)
> a<-table.row.start(a)
> a<-table.element(a, 'R-squared',1,TRUE)
> a<-table.element(a, mysum$r.squared)
> a<-table.row.end(a)
> a<-table.row.start(a)
> a<-table.element(a, 'Adjusted R-squared',1,TRUE)
> a<-table.element(a, mysum$adj.r.squared)
> a<-table.row.end(a)
> a<-table.row.start(a)
> a<-table.element(a, 'F-TEST (value)',1,TRUE)
> a<-table.element(a, mysum$fstatistic[1])
> a<-table.row.end(a)
> a<-table.row.start(a)
> a<-table.element(a, 'F-TEST (DF numerator)',1,TRUE)
> a<-table.element(a, mysum$fstatistic[2])
> a<-table.row.end(a)
> a<-table.row.start(a)
> a<-table.element(a, 'F-TEST (DF denominator)',1,TRUE)
> a<-table.element(a, mysum$fstatistic[3])
> a<-table.row.end(a)
> a<-table.row.start(a)
> a<-table.element(a, 'p-value',1,TRUE)
> a<-table.element(a, 1-pf(mysum$fstatistic[1],mysum$fstatistic[2],mysum$fstatistic[3]))
> a<-table.row.end(a)
> a<-table.row.start(a)
> a<-table.element(a, 'Multiple Linear Regression - Residual Statistics', 2, TRUE)
> a<-table.row.end(a)
> a<-table.row.start(a)
> a<-table.element(a, 'Residual Standard Deviation',1,TRUE)
> a<-table.element(a, mysum$sigma)
> a<-table.row.end(a)
> a<-table.row.start(a)
> a<-table.element(a, 'Sum Squared Residuals',1,TRUE)
> a<-table.element(a, sum(myerror*myerror))
> a<-table.row.end(a)
> a<-table.end(a)
> table.save(a,file="/var/www/html/rcomp/tmp/13744v1258624337.tab")
> a<-table.start()
> a<-table.row.start(a)
> a<-table.element(a, 'Multiple Linear Regression - Actuals, Interpolation, and Residuals', 4, TRUE)
> a<-table.row.end(a)
> a<-table.row.start(a)
> a<-table.element(a, 'Time or Index', 1, TRUE)
> a<-table.element(a, 'Actuals', 1, TRUE)
> a<-table.element(a, 'Interpolation
Forecast', 1, TRUE)
> a<-table.element(a, 'Residuals
Prediction Error', 1, TRUE)
> a<-table.row.end(a)
> for (i in 1:n) {
+ a<-table.row.start(a)
+ a<-table.element(a,i, 1, TRUE)
+ a<-table.element(a,x[i])
+ a<-table.element(a,x[i]-mysum$resid[i])
+ a<-table.element(a,mysum$resid[i])
+ a<-table.row.end(a)
+ }
> a<-table.end(a)
> table.save(a,file="/var/www/html/rcomp/tmp/14u6lj1258624337.tab")
> if (n > n25) {
+ a<-table.start()
+ a<-table.row.start(a)
+ a<-table.element(a,'Goldfeld-Quandt test for Heteroskedasticity',4,TRUE)
+ a<-table.row.end(a)
+ a<-table.row.start(a)
+ a<-table.element(a,'p-values',header=TRUE)
+ a<-table.element(a,'Alternative Hypothesis',3,header=TRUE)
+ a<-table.row.end(a)
+ a<-table.row.start(a)
+ a<-table.element(a,'breakpoint index',header=TRUE)
+ a<-table.element(a,'greater',header=TRUE)
+ a<-table.element(a,'2-sided',header=TRUE)
+ a<-table.element(a,'less',header=TRUE)
+ a<-table.row.end(a)
+ for (mypoint in kp3:nmkm3) {
+ a<-table.row.start(a)
+ a<-table.element(a,mypoint,header=TRUE)
+ a<-table.element(a,gqarr[mypoint-kp3+1,1])
+ a<-table.element(a,gqarr[mypoint-kp3+1,2])
+ a<-table.element(a,gqarr[mypoint-kp3+1,3])
+ a<-table.row.end(a)
+ }
+ a<-table.end(a)
+ table.save(a,file="/var/www/html/rcomp/tmp/15an1v1258624337.tab")
+ a<-table.start()
+ a<-table.row.start(a)
+ a<-table.element(a,'Meta Analysis of Goldfeld-Quandt test for Heteroskedasticity',4,TRUE)
+ a<-table.row.end(a)
+ a<-table.row.start(a)
+ a<-table.element(a,'Description',header=TRUE)
+ a<-table.element(a,'# significant tests',header=TRUE)
+ a<-table.element(a,'% significant tests',header=TRUE)
+ a<-table.element(a,'OK/NOK',header=TRUE)
+ a<-table.row.end(a)
+ a<-table.row.start(a)
+ a<-table.element(a,'1% type I error level',header=TRUE)
+ a<-table.element(a,numsignificant1)
+ a<-table.element(a,numsignificant1/numgqtests)
+ if (numsignificant1/numgqtests < 0.01) dum <- 'OK' else dum <- 'NOK'
+ a<-table.element(a,dum)
+ a<-table.row.end(a)
+ a<-table.row.start(a)
+ a<-table.element(a,'5% type I error level',header=TRUE)
+ a<-table.element(a,numsignificant5)
+ a<-table.element(a,numsignificant5/numgqtests)
+ if (numsignificant5/numgqtests < 0.05) dum <- 'OK' else dum <- 'NOK'
+ a<-table.element(a,dum)
+ a<-table.row.end(a)
+ a<-table.row.start(a)
+ a<-table.element(a,'10% type I error level',header=TRUE)
+ a<-table.element(a,numsignificant10)
+ a<-table.element(a,numsignificant10/numgqtests)
+ if (numsignificant10/numgqtests < 0.1) dum <- 'OK' else dum <- 'NOK'
+ a<-table.element(a,dum)
+ a<-table.row.end(a)
+ a<-table.end(a)
+ table.save(a,file="/var/www/html/rcomp/tmp/16699p1258624337.tab")
+ }
>
> system("convert tmp/1pgqd1258624337.ps tmp/1pgqd1258624337.png")
> system("convert tmp/2h00c1258624337.ps tmp/2h00c1258624337.png")
> system("convert tmp/3fpcv1258624337.ps tmp/3fpcv1258624337.png")
> system("convert tmp/4sodh1258624337.ps tmp/4sodh1258624337.png")
> system("convert tmp/5f7zb1258624337.ps tmp/5f7zb1258624337.png")
> system("convert tmp/64j3g1258624337.ps tmp/64j3g1258624337.png")
> system("convert tmp/7lsm51258624337.ps tmp/7lsm51258624337.png")
> system("convert tmp/8uh5v1258624337.ps tmp/8uh5v1258624337.png")
> system("convert tmp/9395n1258624337.ps tmp/9395n1258624337.png")
> system("convert tmp/10ealw1258624337.ps tmp/10ealw1258624337.png")
>
>
> proc.time()
user system elapsed
2.579 1.595 10.486