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 = '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 t
1 3397 562 1 0 0 0 0 0 0 0 0 0 0 1
2 3971 561 0 1 0 0 0 0 0 0 0 0 0 2
3 4625 555 0 0 1 0 0 0 0 0 0 0 0 3
4 4486 544 0 0 0 1 0 0 0 0 0 0 0 4
5 4132 537 0 0 0 0 1 0 0 0 0 0 0 5
6 4685 543 0 0 0 0 0 1 0 0 0 0 0 6
7 3172 594 0 0 0 0 0 0 1 0 0 0 0 7
8 4280 611 0 0 0 0 0 0 0 1 0 0 0 8
9 4207 613 0 0 0 0 0 0 0 0 1 0 0 9
10 4158 611 0 0 0 0 0 0 0 0 0 1 0 10
11 3933 594 0 0 0 0 0 0 0 0 0 0 1 11
12 3151 595 0 0 0 0 0 0 0 0 0 0 0 12
13 3616 591 1 0 0 0 0 0 0 0 0 0 0 13
14 4221 589 0 1 0 0 0 0 0 0 0 0 0 14
15 4436 584 0 0 1 0 0 0 0 0 0 0 0 15
16 4807 573 0 0 0 1 0 0 0 0 0 0 0 16
17 4849 567 0 0 0 0 1 0 0 0 0 0 0 17
18 5024 569 0 0 0 0 0 1 0 0 0 0 0 18
19 3521 621 0 0 0 0 0 0 1 0 0 0 0 19
20 4650 629 0 0 0 0 0 0 0 1 0 0 0 20
21 5393 628 0 0 0 0 0 0 0 0 1 0 0 21
22 5147 612 0 0 0 0 0 0 0 0 0 1 0 22
23 4845 595 0 0 0 0 0 0 0 0 0 0 1 23
24 3995 597 0 0 0 0 0 0 0 0 0 0 0 24
25 4493 593 1 0 0 0 0 0 0 0 0 0 0 25
26 4680 590 0 1 0 0 0 0 0 0 0 0 0 26
27 5463 580 0 0 1 0 0 0 0 0 0 0 0 27
28 4761 574 0 0 0 1 0 0 0 0 0 0 0 28
29 5307 573 0 0 0 0 1 0 0 0 0 0 0 29
30 5069 573 0 0 0 0 0 1 0 0 0 0 0 30
31 3501 620 0 0 0 0 0 0 1 0 0 0 0 31
32 4952 626 0 0 0 0 0 0 0 1 0 0 0 32
33 5152 620 0 0 0 0 0 0 0 0 1 0 0 33
34 5317 588 0 0 0 0 0 0 0 0 0 1 0 34
35 5189 566 0 0 0 0 0 0 0 0 0 0 1 35
36 4030 557 0 0 0 0 0 0 0 0 0 0 0 36
37 4420 561 1 0 0 0 0 0 0 0 0 0 0 37
38 4571 549 0 1 0 0 0 0 0 0 0 0 0 38
39 4551 532 0 0 1 0 0 0 0 0 0 0 0 39
40 4819 526 0 0 0 1 0 0 0 0 0 0 0 40
41 5133 511 0 0 0 0 1 0 0 0 0 0 0 41
42 4532 499 0 0 0 0 0 1 0 0 0 0 0 42
43 3339 555 0 0 0 0 0 0 1 0 0 0 0 43
44 4380 565 0 0 0 0 0 0 0 1 0 0 0 44
45 4632 542 0 0 0 0 0 0 0 0 1 0 0 45
46 4719 527 0 0 0 0 0 0 0 0 0 1 0 46
47 4212 510 0 0 0 0 0 0 0 0 0 0 1 47
48 3615 514 0 0 0 0 0 0 0 0 0 0 0 48
49 3420 517 1 0 0 0 0 0 0 0 0 0 0 49
50 4571 508 0 1 0 0 0 0 0 0 0 0 0 50
51 4407 493 0 0 1 0 0 0 0 0 0 0 0 51
52 4386 490 0 0 0 1 0 0 0 0 0 0 0 52
53 4386 469 0 0 0 0 1 0 0 0 0 0 0 53
54 4744 478 0 0 0 0 0 1 0 0 0 0 0 54
55 3185 528 0 0 0 0 0 0 1 0 0 0 0 55
56 3890 534 0 0 0 0 0 0 0 1 0 0 0 56
57 4520 518 0 0 0 0 0 0 0 0 1 0 0 57
58 3990 506 0 0 0 0 0 0 0 0 0 1 0 58
59 3809 502 0 0 0 0 0 0 0 0 0 0 1 59
60 3236 516 0 0 0 0 0 0 0 0 0 0 0 60
61 3551 528 1 0 0 0 0 0 0 0 0 0 0 61
62 3264 533 0 1 0 0 0 0 0 0 0 0 0 62
63 3579 536 0 0 1 0 0 0 0 0 0 0 0 63
64 3537 537 0 0 0 1 0 0 0 0 0 0 0 64
65 3038 524 0 0 0 0 1 0 0 0 0 0 0 65
66 2888 536 0 0 0 0 0 1 0 0 0 0 0 66
67 2198 587 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) totWL M1 M2 M3 M4
2727.259 2.095 165.007 577.473 900.048 876.402
M5 M6 M7 M8 M9 M10
914.516 932.698 -504.209 715.265 1092.051 1017.664
M11 t
789.277 -7.951
> (mysum <- summary(mylm))
Call:
lm(formula = df)
Residuals:
Min 1Q Median 3Q Max
-1370.1 -292.8 130.1 305.7 835.3
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 2727.259 1627.252 1.676 0.09963 .
totWL 2.095 2.662 0.787 0.43485
M1 165.007 329.987 0.500 0.61912
M2 577.473 330.069 1.750 0.08598 .
M3 900.048 331.377 2.716 0.00890 **
M4 876.402 332.831 2.633 0.01106 *
M5 914.516 337.125 2.713 0.00898 **
M6 932.698 334.940 2.785 0.00741 **
M7 -504.209 338.747 -1.488 0.14256
M8 715.265 354.705 2.017 0.04883 *
M9 1092.051 350.325 3.117 0.00295 **
M10 1017.664 345.223 2.948 0.00475 **
M11 789.277 344.199 2.293 0.02584 *
t -7.951 4.885 -1.628 0.10955
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 544 on 53 degrees of freedom
Multiple R-squared: 0.5249, Adjusted R-squared: 0.4084
F-statistic: 4.504 on 13 and 53 DF, p-value: 4.166e-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.31527503 0.6305501 0.6847250
[2,] 0.26982649 0.5396530 0.7301735
[3,] 0.19189008 0.3837802 0.8081099
[4,] 0.13241690 0.2648338 0.8675831
[5,] 0.21691017 0.4338203 0.7830898
[6,] 0.17438928 0.3487786 0.8256107
[7,] 0.12996839 0.2599368 0.8700316
[8,] 0.10378427 0.2075685 0.8962157
[9,] 0.06972687 0.1394537 0.9302731
[10,] 0.13046944 0.2609389 0.8695306
[11,] 0.09778312 0.1955662 0.9022169
[12,] 0.37138660 0.7427732 0.6286134
[13,] 0.29536391 0.5907278 0.7046361
[14,] 0.35980477 0.7196095 0.6401952
[15,] 0.48060782 0.9612156 0.5193922
[16,] 0.43307922 0.8661584 0.5669208
[17,] 0.40260945 0.8052189 0.5973905
[18,] 0.37147608 0.7429522 0.6285239
[19,] 0.44653621 0.8930724 0.5534638
[20,] 0.38796682 0.7759336 0.6120332
[21,] 0.37992849 0.7598570 0.6200715
[22,] 0.36873509 0.7374702 0.6312649
[23,] 0.50016666 0.9996667 0.4998333
[24,] 0.42845545 0.8569109 0.5715446
[25,] 0.59452136 0.8109573 0.4054786
[26,] 0.58129205 0.8374159 0.4187079
[27,] 0.50505167 0.9898967 0.4949483
[28,] 0.46753843 0.9350769 0.5324616
[29,] 0.36574229 0.7314846 0.6342577
[30,] 0.41550108 0.8310022 0.5844989
[31,] 0.35465077 0.7093015 0.6453492
[32,] 0.26311968 0.5262394 0.7368803
[33,] 0.42621206 0.8524241 0.5737879
[34,] 0.51124950 0.9775010 0.4887505
> postscript(file="/var/www/html/rcomp/tmp/1vug61261222471.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/27dll1261222471.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/34pu51261222471.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/4i9on1261222471.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/56zp01261222472.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
-664.6756503 -493.0955094 -141.1503861 -225.5089836 -595.0079653
6 7 8 9 10
-64.8086314 -239.7928068 -378.9297101 -824.9551395 -787.4274315
11 12 13 14 15
-740.4755160 -727.3425900 -411.0187875 -206.3436994 -295.4935233
16 17 18 19 20
130.1478792 154.5539503 315.1330729 148.0539503 48.7715715
21 22 23 24 25
425.0309836 294.8879519 264.8398674 207.8778462 557.2016487
26 27 28 29 30
345.9716839 835.2965959 177.4632626 695.3945978 447.1636148
31 32 33 34 35
225.5592280 452.4667435 296.2008915 610.5770144 765.0036657
36 37 38 39 40
422.0860634 646.6502886 418.2748483 119.2643904 431.4310570
41 42 43 44 45
746.6916526 160.6000356 295.1411243 103.6688512 35.0171009
46 47 48 49 50
235.7791221 0.7310375 192.5791221 -165.7617056 599.5780126
51 52 53 54 55
152.3776604 169.2594855 183.0897641 512.0042566 293.1150283
56 57 58 59 60
-225.9774562 68.7061634 -353.8166569 -290.0990546 -95.2004417
61 62 63 64 65
37.6042061 -664.3853360 -670.2947372 -682.7927007 -1184.7219994
66 67
-1370.0923485 -722.0765239
> postscript(file="/var/www/html/rcomp/tmp/6cqe11261222472.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 -664.6756503 NA
1 -493.0955094 -664.6756503
2 -141.1503861 -493.0955094
3 -225.5089836 -141.1503861
4 -595.0079653 -225.5089836
5 -64.8086314 -595.0079653
6 -239.7928068 -64.8086314
7 -378.9297101 -239.7928068
8 -824.9551395 -378.9297101
9 -787.4274315 -824.9551395
10 -740.4755160 -787.4274315
11 -727.3425900 -740.4755160
12 -411.0187875 -727.3425900
13 -206.3436994 -411.0187875
14 -295.4935233 -206.3436994
15 130.1478792 -295.4935233
16 154.5539503 130.1478792
17 315.1330729 154.5539503
18 148.0539503 315.1330729
19 48.7715715 148.0539503
20 425.0309836 48.7715715
21 294.8879519 425.0309836
22 264.8398674 294.8879519
23 207.8778462 264.8398674
24 557.2016487 207.8778462
25 345.9716839 557.2016487
26 835.2965959 345.9716839
27 177.4632626 835.2965959
28 695.3945978 177.4632626
29 447.1636148 695.3945978
30 225.5592280 447.1636148
31 452.4667435 225.5592280
32 296.2008915 452.4667435
33 610.5770144 296.2008915
34 765.0036657 610.5770144
35 422.0860634 765.0036657
36 646.6502886 422.0860634
37 418.2748483 646.6502886
38 119.2643904 418.2748483
39 431.4310570 119.2643904
40 746.6916526 431.4310570
41 160.6000356 746.6916526
42 295.1411243 160.6000356
43 103.6688512 295.1411243
44 35.0171009 103.6688512
45 235.7791221 35.0171009
46 0.7310375 235.7791221
47 192.5791221 0.7310375
48 -165.7617056 192.5791221
49 599.5780126 -165.7617056
50 152.3776604 599.5780126
51 169.2594855 152.3776604
52 183.0897641 169.2594855
53 512.0042566 183.0897641
54 293.1150283 512.0042566
55 -225.9774562 293.1150283
56 68.7061634 -225.9774562
57 -353.8166569 68.7061634
58 -290.0990546 -353.8166569
59 -95.2004417 -290.0990546
60 37.6042061 -95.2004417
61 -664.3853360 37.6042061
62 -670.2947372 -664.3853360
63 -682.7927007 -670.2947372
64 -1184.7219994 -682.7927007
65 -1370.0923485 -1184.7219994
66 -722.0765239 -1370.0923485
67 NA -722.0765239
> dum1 <- dum[2:length(myerror),]
> dum1
lag(myerror, k = 1) myerror
[1,] -493.0955094 -664.6756503
[2,] -141.1503861 -493.0955094
[3,] -225.5089836 -141.1503861
[4,] -595.0079653 -225.5089836
[5,] -64.8086314 -595.0079653
[6,] -239.7928068 -64.8086314
[7,] -378.9297101 -239.7928068
[8,] -824.9551395 -378.9297101
[9,] -787.4274315 -824.9551395
[10,] -740.4755160 -787.4274315
[11,] -727.3425900 -740.4755160
[12,] -411.0187875 -727.3425900
[13,] -206.3436994 -411.0187875
[14,] -295.4935233 -206.3436994
[15,] 130.1478792 -295.4935233
[16,] 154.5539503 130.1478792
[17,] 315.1330729 154.5539503
[18,] 148.0539503 315.1330729
[19,] 48.7715715 148.0539503
[20,] 425.0309836 48.7715715
[21,] 294.8879519 425.0309836
[22,] 264.8398674 294.8879519
[23,] 207.8778462 264.8398674
[24,] 557.2016487 207.8778462
[25,] 345.9716839 557.2016487
[26,] 835.2965959 345.9716839
[27,] 177.4632626 835.2965959
[28,] 695.3945978 177.4632626
[29,] 447.1636148 695.3945978
[30,] 225.5592280 447.1636148
[31,] 452.4667435 225.5592280
[32,] 296.2008915 452.4667435
[33,] 610.5770144 296.2008915
[34,] 765.0036657 610.5770144
[35,] 422.0860634 765.0036657
[36,] 646.6502886 422.0860634
[37,] 418.2748483 646.6502886
[38,] 119.2643904 418.2748483
[39,] 431.4310570 119.2643904
[40,] 746.6916526 431.4310570
[41,] 160.6000356 746.6916526
[42,] 295.1411243 160.6000356
[43,] 103.6688512 295.1411243
[44,] 35.0171009 103.6688512
[45,] 235.7791221 35.0171009
[46,] 0.7310375 235.7791221
[47,] 192.5791221 0.7310375
[48,] -165.7617056 192.5791221
[49,] 599.5780126 -165.7617056
[50,] 152.3776604 599.5780126
[51,] 169.2594855 152.3776604
[52,] 183.0897641 169.2594855
[53,] 512.0042566 183.0897641
[54,] 293.1150283 512.0042566
[55,] -225.9774562 293.1150283
[56,] 68.7061634 -225.9774562
[57,] -353.8166569 68.7061634
[58,] -290.0990546 -353.8166569
[59,] -95.2004417 -290.0990546
[60,] 37.6042061 -95.2004417
[61,] -664.3853360 37.6042061
[62,] -670.2947372 -664.3853360
[63,] -682.7927007 -670.2947372
[64,] -1184.7219994 -682.7927007
[65,] -1370.0923485 -1184.7219994
[66,] -722.0765239 -1370.0923485
> z <- as.data.frame(dum1)
> z
lag(myerror, k = 1) myerror
1 -493.0955094 -664.6756503
2 -141.1503861 -493.0955094
3 -225.5089836 -141.1503861
4 -595.0079653 -225.5089836
5 -64.8086314 -595.0079653
6 -239.7928068 -64.8086314
7 -378.9297101 -239.7928068
8 -824.9551395 -378.9297101
9 -787.4274315 -824.9551395
10 -740.4755160 -787.4274315
11 -727.3425900 -740.4755160
12 -411.0187875 -727.3425900
13 -206.3436994 -411.0187875
14 -295.4935233 -206.3436994
15 130.1478792 -295.4935233
16 154.5539503 130.1478792
17 315.1330729 154.5539503
18 148.0539503 315.1330729
19 48.7715715 148.0539503
20 425.0309836 48.7715715
21 294.8879519 425.0309836
22 264.8398674 294.8879519
23 207.8778462 264.8398674
24 557.2016487 207.8778462
25 345.9716839 557.2016487
26 835.2965959 345.9716839
27 177.4632626 835.2965959
28 695.3945978 177.4632626
29 447.1636148 695.3945978
30 225.5592280 447.1636148
31 452.4667435 225.5592280
32 296.2008915 452.4667435
33 610.5770144 296.2008915
34 765.0036657 610.5770144
35 422.0860634 765.0036657
36 646.6502886 422.0860634
37 418.2748483 646.6502886
38 119.2643904 418.2748483
39 431.4310570 119.2643904
40 746.6916526 431.4310570
41 160.6000356 746.6916526
42 295.1411243 160.6000356
43 103.6688512 295.1411243
44 35.0171009 103.6688512
45 235.7791221 35.0171009
46 0.7310375 235.7791221
47 192.5791221 0.7310375
48 -165.7617056 192.5791221
49 599.5780126 -165.7617056
50 152.3776604 599.5780126
51 169.2594855 152.3776604
52 183.0897641 169.2594855
53 512.0042566 183.0897641
54 293.1150283 512.0042566
55 -225.9774562 293.1150283
56 68.7061634 -225.9774562
57 -353.8166569 68.7061634
58 -290.0990546 -353.8166569
59 -95.2004417 -290.0990546
60 37.6042061 -95.2004417
61 -664.3853360 37.6042061
62 -670.2947372 -664.3853360
63 -682.7927007 -670.2947372
64 -1184.7219994 -682.7927007
65 -1370.0923485 -1184.7219994
66 -722.0765239 -1370.0923485
> 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/7sd561261222472.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/8ul7x1261222472.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/9fh8w1261222472.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/10ti9r1261222472.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/11qw8t1261222472.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/12zqso1261222472.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/137cwo1261222472.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/1407wd1261222472.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/15d87l1261222472.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/160om11261222472.tab")
+ }
> try(system("convert tmp/1vug61261222471.ps tmp/1vug61261222471.png",intern=TRUE))
character(0)
> try(system("convert tmp/27dll1261222471.ps tmp/27dll1261222471.png",intern=TRUE))
character(0)
> try(system("convert tmp/34pu51261222471.ps tmp/34pu51261222471.png",intern=TRUE))
character(0)
> try(system("convert tmp/4i9on1261222471.ps tmp/4i9on1261222471.png",intern=TRUE))
character(0)
> try(system("convert tmp/56zp01261222472.ps tmp/56zp01261222472.png",intern=TRUE))
character(0)
> try(system("convert tmp/6cqe11261222472.ps tmp/6cqe11261222472.png",intern=TRUE))
character(0)
> try(system("convert tmp/7sd561261222472.ps tmp/7sd561261222472.png",intern=TRUE))
character(0)
> try(system("convert tmp/8ul7x1261222472.ps tmp/8ul7x1261222472.png",intern=TRUE))
character(0)
> try(system("convert tmp/9fh8w1261222472.ps tmp/9fh8w1261222472.png",intern=TRUE))
character(0)
> try(system("convert tmp/10ti9r1261222472.ps tmp/10ti9r1261222472.png",intern=TRUE))
character(0)
>
>
> proc.time()
user system elapsed
2.472 1.549 3.421