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(3397,562,3971,561,4625,555,4486,544,4132,537,4685,543,3172,594,4280,611,4207,613,4158,611,3933,594,3151,595,3616,591,4221,589,4436,584,4807,573,4849,567,5024,569,3521,621,4650,629,5393,628,5147,612,4845,595,3995,597,4493,593,4680,590,5463,580,4761,574,5307,573,5069,573,3501,620,4952,626,5152,620,5317,588,5189,566,4030,557,4420,561,4571,549,4551,532,4819,526,5133,511,4532,499,3339,555,4380,565,4632,542,4719,527,4212,510,3615,514,3420,517,4571,508,4407,493,4386,490,4386,469,4744,478,3185,528,3890,534,4520,518,3990,506,3809,502,3236,516,3551,528,3264,533,3579,536,3537,537,3038,524,2888,536,2198,587),dim=c(2,67),dimnames=list(c('wng','totWL'),1:67))
> y <- array(NA,dim=c(2,67),dimnames=list(c('wng','totWL'),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 = 'No 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
wng totWL M1 M2 M3 M4 M5 M6 M7 M8 M9 M10 M11
1 3397 562 1 0 0 0 0 0 0 0 0 0 0
2 3971 561 0 1 0 0 0 0 0 0 0 0 0
3 4625 555 0 0 1 0 0 0 0 0 0 0 0
4 4486 544 0 0 0 1 0 0 0 0 0 0 0
5 4132 537 0 0 0 0 1 0 0 0 0 0 0
6 4685 543 0 0 0 0 0 1 0 0 0 0 0
7 3172 594 0 0 0 0 0 0 1 0 0 0 0
8 4280 611 0 0 0 0 0 0 0 1 0 0 0
9 4207 613 0 0 0 0 0 0 0 0 1 0 0
10 4158 611 0 0 0 0 0 0 0 0 0 1 0
11 3933 594 0 0 0 0 0 0 0 0 0 0 1
12 3151 595 0 0 0 0 0 0 0 0 0 0 0
13 3616 591 1 0 0 0 0 0 0 0 0 0 0
14 4221 589 0 1 0 0 0 0 0 0 0 0 0
15 4436 584 0 0 1 0 0 0 0 0 0 0 0
16 4807 573 0 0 0 1 0 0 0 0 0 0 0
17 4849 567 0 0 0 0 1 0 0 0 0 0 0
18 5024 569 0 0 0 0 0 1 0 0 0 0 0
19 3521 621 0 0 0 0 0 0 1 0 0 0 0
20 4650 629 0 0 0 0 0 0 0 1 0 0 0
21 5393 628 0 0 0 0 0 0 0 0 1 0 0
22 5147 612 0 0 0 0 0 0 0 0 0 1 0
23 4845 595 0 0 0 0 0 0 0 0 0 0 1
24 3995 597 0 0 0 0 0 0 0 0 0 0 0
25 4493 593 1 0 0 0 0 0 0 0 0 0 0
26 4680 590 0 1 0 0 0 0 0 0 0 0 0
27 5463 580 0 0 1 0 0 0 0 0 0 0 0
28 4761 574 0 0 0 1 0 0 0 0 0 0 0
29 5307 573 0 0 0 0 1 0 0 0 0 0 0
30 5069 573 0 0 0 0 0 1 0 0 0 0 0
31 3501 620 0 0 0 0 0 0 1 0 0 0 0
32 4952 626 0 0 0 0 0 0 0 1 0 0 0
33 5152 620 0 0 0 0 0 0 0 0 1 0 0
34 5317 588 0 0 0 0 0 0 0 0 0 1 0
35 5189 566 0 0 0 0 0 0 0 0 0 0 1
36 4030 557 0 0 0 0 0 0 0 0 0 0 0
37 4420 561 1 0 0 0 0 0 0 0 0 0 0
38 4571 549 0 1 0 0 0 0 0 0 0 0 0
39 4551 532 0 0 1 0 0 0 0 0 0 0 0
40 4819 526 0 0 0 1 0 0 0 0 0 0 0
41 5133 511 0 0 0 0 1 0 0 0 0 0 0
42 4532 499 0 0 0 0 0 1 0 0 0 0 0
43 3339 555 0 0 0 0 0 0 1 0 0 0 0
44 4380 565 0 0 0 0 0 0 0 1 0 0 0
45 4632 542 0 0 0 0 0 0 0 0 1 0 0
46 4719 527 0 0 0 0 0 0 0 0 0 1 0
47 4212 510 0 0 0 0 0 0 0 0 0 0 1
48 3615 514 0 0 0 0 0 0 0 0 0 0 0
49 3420 517 1 0 0 0 0 0 0 0 0 0 0
50 4571 508 0 1 0 0 0 0 0 0 0 0 0
51 4407 493 0 0 1 0 0 0 0 0 0 0 0
52 4386 490 0 0 0 1 0 0 0 0 0 0 0
53 4386 469 0 0 0 0 1 0 0 0 0 0 0
54 4744 478 0 0 0 0 0 1 0 0 0 0 0
55 3185 528 0 0 0 0 0 0 1 0 0 0 0
56 3890 534 0 0 0 0 0 0 0 1 0 0 0
57 4520 518 0 0 0 0 0 0 0 0 1 0 0
58 3990 506 0 0 0 0 0 0 0 0 0 1 0
59 3809 502 0 0 0 0 0 0 0 0 0 0 1
60 3236 516 0 0 0 0 0 0 0 0 0 0 0
61 3551 528 1 0 0 0 0 0 0 0 0 0 0
62 3264 533 0 1 0 0 0 0 0 0 0 0 0
63 3579 536 0 0 1 0 0 0 0 0 0 0 0
64 3537 537 0 0 0 1 0 0 0 0 0 0 0
65 3038 524 0 0 0 0 1 0 0 0 0 0 0
66 2888 536 0 0 0 0 0 1 0 0 0 0 0
67 2198 587 0 0 0 0 0 0 1 0 0 0 0
> k <- length(x[1,])
> df <- as.data.frame(x)
> (mylm <- lm(df))
Call:
lm(formula = df)
Coefficients:
(Intercept) totWL M1 M2 M3 M4
736.698 5.161 195.971 611.729 951.907 938.709
M5 M6 M7 M8 M9 M10
1001.070 1002.613 -599.145 632.996 1028.816 993.702
M11
804.587
> (mysum <- summary(mylm))
Call:
lm(formula = df)
Residuals:
Min 1Q Median 3Q Max
-1617.8 -255.1 143.1 329.6 780.8
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 736.698 1089.649 0.676 0.50187
totWL 5.161 1.909 2.703 0.00917 **
M1 195.971 334.430 0.586 0.56033
M2 611.729 334.389 1.829 0.07286 .
M3 951.907 334.840 2.843 0.00630 **
M4 938.709 335.632 2.797 0.00714 **
M5 1001.070 337.949 2.962 0.00453 **
M6 1002.613 337.208 2.973 0.00440 **
M7 -599.145 338.744 -1.769 0.08259 .
M8 632.996 356.405 1.776 0.08136 .
M9 1028.816 353.439 2.911 0.00523 **
M10 993.702 350.135 2.838 0.00638 **
M11 804.587 349.284 2.304 0.02513 *
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 552.2 on 54 degrees of freedom
Multiple R-squared: 0.5011, Adjusted R-squared: 0.3903
F-statistic: 4.521 on 12 and 54 DF, p-value: 5.37e-05
> 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.029816635 0.059633270 0.9701834
[2,] 0.038105915 0.076211831 0.9618941
[3,] 0.012922098 0.025844196 0.9870779
[4,] 0.003963676 0.007927352 0.9960363
[5,] 0.001383811 0.002767623 0.9986162
[6,] 0.020896502 0.041793003 0.9791035
[7,] 0.048160628 0.096321256 0.9518394
[8,] 0.064841122 0.129682243 0.9351589
[9,] 0.068802758 0.137605516 0.9311972
[10,] 0.073402150 0.146804300 0.9265979
[11,] 0.048867908 0.097735816 0.9511321
[12,] 0.059738720 0.119477440 0.9402613
[13,] 0.036795670 0.073591340 0.9632043
[14,] 0.029283589 0.058567177 0.9707164
[15,] 0.019644384 0.039288768 0.9803556
[16,] 0.011358843 0.022717686 0.9886412
[17,] 0.008063186 0.016126371 0.9919368
[18,] 0.005294717 0.010589434 0.9947053
[19,] 0.014804759 0.029609518 0.9851952
[20,] 0.054756316 0.109512632 0.9452437
[21,] 0.071006552 0.142013104 0.9289934
[22,] 0.124163637 0.248327275 0.8758364
[23,] 0.143521245 0.287042490 0.8564788
[24,] 0.133488071 0.266976141 0.8665119
[25,] 0.161904102 0.323808204 0.8380959
[26,] 0.481587884 0.963175768 0.5184121
[27,] 0.454580511 0.909161022 0.5454195
[28,] 0.434095590 0.868191179 0.5659044
[29,] 0.561033034 0.877933932 0.4389670
[30,] 0.539168922 0.921662156 0.4608311
[31,] 0.808642324 0.382715353 0.1913577
[32,] 0.822034839 0.355930321 0.1779652
[33,] 0.761287939 0.477424121 0.2387121
[34,] 0.709764808 0.580470384 0.2902352
[35,] 0.837039903 0.325920195 0.1629601
[36,] 0.708989569 0.582020862 0.2910104
> postscript(file="/var/www/html/rcomp/tmp/1dwnc1261154142.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/2jo761261154142.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/32cd11261154142.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/4zs3x1261154142.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/5qywl1261154142.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 = 67
Frequency = 1
1 2 3 4 5 6
-436.371308 -272.968355 71.821730 2.795359 -377.436182 143.052742
7 8 9 10 11 12
-31.420359 -243.305064 -722.448102 -726.010761 -674.152533 -656.726584
13 14 15 16 17 18
-367.051689 -167.487343 -266.858651 174.114978 184.722045 347.856539
19 20 21 22 23 24
178.222045 33.789872 386.131011 257.827846 232.686074 176.950631
25 26 27 28 29 30
499.625526 286.351264 780.786919 122.953585 611.753690 372.210969
31 32 33 34 35 36
163.383437 351.274049 186.422151 551.701265 726.366455 418.406329
37 38 39 40 41 42
591.790084 388.968355 116.533756 428.700422 757.760022 217.154010
43 44 45 46 47 48
336.873946 94.118988 69.010761 268.546204 38.404432 225.346204
49 50 51 52 53 54
-181.108648 600.585445 173.828061 181.510550 227.538504 537.543251
55 56 57 58 59 60
322.231542 -235.877846 80.884180 -352.064555 -323.304428 -163.976581
61 62 63 64 65 66
-106.883965 -835.449366 -876.111814 -910.074894 -1404.338080 -1617.817511
67
-969.290612
> postscript(file="/var/www/html/rcomp/tmp/64pvz1261154142.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 = 67
Frequency = 1
lag(myerror, k = 1) myerror
0 -436.371308 NA
1 -272.968355 -436.371308
2 71.821730 -272.968355
3 2.795359 71.821730
4 -377.436182 2.795359
5 143.052742 -377.436182
6 -31.420359 143.052742
7 -243.305064 -31.420359
8 -722.448102 -243.305064
9 -726.010761 -722.448102
10 -674.152533 -726.010761
11 -656.726584 -674.152533
12 -367.051689 -656.726584
13 -167.487343 -367.051689
14 -266.858651 -167.487343
15 174.114978 -266.858651
16 184.722045 174.114978
17 347.856539 184.722045
18 178.222045 347.856539
19 33.789872 178.222045
20 386.131011 33.789872
21 257.827846 386.131011
22 232.686074 257.827846
23 176.950631 232.686074
24 499.625526 176.950631
25 286.351264 499.625526
26 780.786919 286.351264
27 122.953585 780.786919
28 611.753690 122.953585
29 372.210969 611.753690
30 163.383437 372.210969
31 351.274049 163.383437
32 186.422151 351.274049
33 551.701265 186.422151
34 726.366455 551.701265
35 418.406329 726.366455
36 591.790084 418.406329
37 388.968355 591.790084
38 116.533756 388.968355
39 428.700422 116.533756
40 757.760022 428.700422
41 217.154010 757.760022
42 336.873946 217.154010
43 94.118988 336.873946
44 69.010761 94.118988
45 268.546204 69.010761
46 38.404432 268.546204
47 225.346204 38.404432
48 -181.108648 225.346204
49 600.585445 -181.108648
50 173.828061 600.585445
51 181.510550 173.828061
52 227.538504 181.510550
53 537.543251 227.538504
54 322.231542 537.543251
55 -235.877846 322.231542
56 80.884180 -235.877846
57 -352.064555 80.884180
58 -323.304428 -352.064555
59 -163.976581 -323.304428
60 -106.883965 -163.976581
61 -835.449366 -106.883965
62 -876.111814 -835.449366
63 -910.074894 -876.111814
64 -1404.338080 -910.074894
65 -1617.817511 -1404.338080
66 -969.290612 -1617.817511
67 NA -969.290612
> dum1 <- dum[2:length(myerror),]
> dum1
lag(myerror, k = 1) myerror
[1,] -272.968355 -436.371308
[2,] 71.821730 -272.968355
[3,] 2.795359 71.821730
[4,] -377.436182 2.795359
[5,] 143.052742 -377.436182
[6,] -31.420359 143.052742
[7,] -243.305064 -31.420359
[8,] -722.448102 -243.305064
[9,] -726.010761 -722.448102
[10,] -674.152533 -726.010761
[11,] -656.726584 -674.152533
[12,] -367.051689 -656.726584
[13,] -167.487343 -367.051689
[14,] -266.858651 -167.487343
[15,] 174.114978 -266.858651
[16,] 184.722045 174.114978
[17,] 347.856539 184.722045
[18,] 178.222045 347.856539
[19,] 33.789872 178.222045
[20,] 386.131011 33.789872
[21,] 257.827846 386.131011
[22,] 232.686074 257.827846
[23,] 176.950631 232.686074
[24,] 499.625526 176.950631
[25,] 286.351264 499.625526
[26,] 780.786919 286.351264
[27,] 122.953585 780.786919
[28,] 611.753690 122.953585
[29,] 372.210969 611.753690
[30,] 163.383437 372.210969
[31,] 351.274049 163.383437
[32,] 186.422151 351.274049
[33,] 551.701265 186.422151
[34,] 726.366455 551.701265
[35,] 418.406329 726.366455
[36,] 591.790084 418.406329
[37,] 388.968355 591.790084
[38,] 116.533756 388.968355
[39,] 428.700422 116.533756
[40,] 757.760022 428.700422
[41,] 217.154010 757.760022
[42,] 336.873946 217.154010
[43,] 94.118988 336.873946
[44,] 69.010761 94.118988
[45,] 268.546204 69.010761
[46,] 38.404432 268.546204
[47,] 225.346204 38.404432
[48,] -181.108648 225.346204
[49,] 600.585445 -181.108648
[50,] 173.828061 600.585445
[51,] 181.510550 173.828061
[52,] 227.538504 181.510550
[53,] 537.543251 227.538504
[54,] 322.231542 537.543251
[55,] -235.877846 322.231542
[56,] 80.884180 -235.877846
[57,] -352.064555 80.884180
[58,] -323.304428 -352.064555
[59,] -163.976581 -323.304428
[60,] -106.883965 -163.976581
[61,] -835.449366 -106.883965
[62,] -876.111814 -835.449366
[63,] -910.074894 -876.111814
[64,] -1404.338080 -910.074894
[65,] -1617.817511 -1404.338080
[66,] -969.290612 -1617.817511
> z <- as.data.frame(dum1)
> z
lag(myerror, k = 1) myerror
1 -272.968355 -436.371308
2 71.821730 -272.968355
3 2.795359 71.821730
4 -377.436182 2.795359
5 143.052742 -377.436182
6 -31.420359 143.052742
7 -243.305064 -31.420359
8 -722.448102 -243.305064
9 -726.010761 -722.448102
10 -674.152533 -726.010761
11 -656.726584 -674.152533
12 -367.051689 -656.726584
13 -167.487343 -367.051689
14 -266.858651 -167.487343
15 174.114978 -266.858651
16 184.722045 174.114978
17 347.856539 184.722045
18 178.222045 347.856539
19 33.789872 178.222045
20 386.131011 33.789872
21 257.827846 386.131011
22 232.686074 257.827846
23 176.950631 232.686074
24 499.625526 176.950631
25 286.351264 499.625526
26 780.786919 286.351264
27 122.953585 780.786919
28 611.753690 122.953585
29 372.210969 611.753690
30 163.383437 372.210969
31 351.274049 163.383437
32 186.422151 351.274049
33 551.701265 186.422151
34 726.366455 551.701265
35 418.406329 726.366455
36 591.790084 418.406329
37 388.968355 591.790084
38 116.533756 388.968355
39 428.700422 116.533756
40 757.760022 428.700422
41 217.154010 757.760022
42 336.873946 217.154010
43 94.118988 336.873946
44 69.010761 94.118988
45 268.546204 69.010761
46 38.404432 268.546204
47 225.346204 38.404432
48 -181.108648 225.346204
49 600.585445 -181.108648
50 173.828061 600.585445
51 181.510550 173.828061
52 227.538504 181.510550
53 537.543251 227.538504
54 322.231542 537.543251
55 -235.877846 322.231542
56 80.884180 -235.877846
57 -352.064555 80.884180
58 -323.304428 -352.064555
59 -163.976581 -323.304428
60 -106.883965 -163.976581
61 -835.449366 -106.883965
62 -876.111814 -835.449366
63 -910.074894 -876.111814
64 -1404.338080 -910.074894
65 -1617.817511 -1404.338080
66 -969.290612 -1617.817511
> 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/7iy471261154142.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/84xq01261154142.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/9ood41261154142.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/10a3sv1261154142.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/11m24z1261154142.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/12gub51261154142.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/13phjr1261154142.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/14220y1261154142.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/15otx61261154142.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/169agz1261154143.tab")
+ }
> try(system("convert tmp/1dwnc1261154142.ps tmp/1dwnc1261154142.png",intern=TRUE))
character(0)
> try(system("convert tmp/2jo761261154142.ps tmp/2jo761261154142.png",intern=TRUE))
character(0)
> try(system("convert tmp/32cd11261154142.ps tmp/32cd11261154142.png",intern=TRUE))
character(0)
> try(system("convert tmp/4zs3x1261154142.ps tmp/4zs3x1261154142.png",intern=TRUE))
character(0)
> try(system("convert tmp/5qywl1261154142.ps tmp/5qywl1261154142.png",intern=TRUE))
character(0)
> try(system("convert tmp/64pvz1261154142.ps tmp/64pvz1261154142.png",intern=TRUE))
character(0)
> try(system("convert tmp/7iy471261154142.ps tmp/7iy471261154142.png",intern=TRUE))
character(0)
> try(system("convert tmp/84xq01261154142.ps tmp/84xq01261154142.png",intern=TRUE))
character(0)
> try(system("convert tmp/9ood41261154142.ps tmp/9ood41261154142.png",intern=TRUE))
character(0)
> try(system("convert tmp/10a3sv1261154142.ps tmp/10a3sv1261154142.png",intern=TRUE))
character(0)
>
>
> proc.time()
user system elapsed
2.451 1.543 3.114