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(452
+ ,1870
+ ,449
+ ,441
+ ,462
+ ,2263
+ ,452
+ ,449
+ ,455
+ ,1802
+ ,462
+ ,452
+ ,461
+ ,1863
+ ,455
+ ,462
+ ,461
+ ,1989
+ ,461
+ ,455
+ ,463
+ ,2197
+ ,461
+ ,461
+ ,462
+ ,2409
+ ,463
+ ,461
+ ,456
+ ,2502
+ ,462
+ ,463
+ ,455
+ ,2593
+ ,456
+ ,462
+ ,456
+ ,2598
+ ,455
+ ,456
+ ,472
+ ,2053
+ ,456
+ ,455
+ ,472
+ ,2213
+ ,472
+ ,456
+ ,471
+ ,2238
+ ,472
+ ,472
+ ,465
+ ,2359
+ ,471
+ ,472
+ ,459
+ ,2151
+ ,465
+ ,471
+ ,465
+ ,2474
+ ,459
+ ,465
+ ,468
+ ,3079
+ ,465
+ ,459
+ ,467
+ ,2312
+ ,468
+ ,465
+ ,463
+ ,2565
+ ,467
+ ,468
+ ,460
+ ,1972
+ ,463
+ ,467
+ ,462
+ ,2484
+ ,460
+ ,463
+ ,461
+ ,2202
+ ,462
+ ,460
+ ,476
+ ,2151
+ ,461
+ ,462
+ ,476
+ ,1976
+ ,476
+ ,461
+ ,471
+ ,2012
+ ,476
+ ,476
+ ,453
+ ,2114
+ ,471
+ ,476
+ ,443
+ ,1772
+ ,453
+ ,471
+ ,442
+ ,1957
+ ,443
+ ,453
+ ,444
+ ,2070
+ ,442
+ ,443
+ ,438
+ ,1990
+ ,444
+ ,442
+ ,427
+ ,2182
+ ,438
+ ,444
+ ,424
+ ,2008
+ ,427
+ ,438
+ ,416
+ ,1916
+ ,424
+ ,427
+ ,406
+ ,2397
+ ,416
+ ,424
+ ,431
+ ,2114
+ ,406
+ ,416
+ ,434
+ ,1778
+ ,431
+ ,406
+ ,418
+ ,1641
+ ,434
+ ,431
+ ,412
+ ,2186
+ ,418
+ ,434
+ ,404
+ ,1773
+ ,412
+ ,418
+ ,409
+ ,1785
+ ,404
+ ,412
+ ,412
+ ,2217
+ ,409
+ ,404
+ ,406
+ ,2153
+ ,412
+ ,409
+ ,398
+ ,1895
+ ,406
+ ,412
+ ,397
+ ,2475
+ ,398
+ ,406
+ ,385
+ ,1793
+ ,397
+ ,398
+ ,390
+ ,2308
+ ,385
+ ,397
+ ,413
+ ,2051
+ ,390
+ ,385
+ ,413
+ ,1898
+ ,413
+ ,390
+ ,401
+ ,2142
+ ,413
+ ,413
+ ,397
+ ,1874
+ ,401
+ ,413
+ ,397
+ ,1560
+ ,397
+ ,401
+ ,409
+ ,1808
+ ,397
+ ,397
+ ,419
+ ,1575
+ ,409
+ ,397
+ ,424
+ ,1525
+ ,419
+ ,409
+ ,428
+ ,1997
+ ,424
+ ,419
+ ,430
+ ,1753
+ ,428
+ ,424
+ ,424
+ ,1623
+ ,430
+ ,428
+ ,433
+ ,2251
+ ,424
+ ,430
+ ,456
+ ,1890
+ ,433
+ ,424)
+ ,dim=c(4
+ ,59)
+ ,dimnames=list(c('wkl'
+ ,'bvg'
+ ,'Y1'
+ ,'Y2')
+ ,1:59))
> y <- array(NA,dim=c(4,59),dimnames=list(c('wkl','bvg','Y1','Y2'),1:59))
> 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
wkl bvg Y1 Y2 M1 M2 M3 M4 M5 M6 M7 M8 M9 M10 M11 t
1 452 1870 449 441 1 0 0 0 0 0 0 0 0 0 0 1
2 462 2263 452 449 0 1 0 0 0 0 0 0 0 0 0 2
3 455 1802 462 452 0 0 1 0 0 0 0 0 0 0 0 3
4 461 1863 455 462 0 0 0 1 0 0 0 0 0 0 0 4
5 461 1989 461 455 0 0 0 0 1 0 0 0 0 0 0 5
6 463 2197 461 461 0 0 0 0 0 1 0 0 0 0 0 6
7 462 2409 463 461 0 0 0 0 0 0 1 0 0 0 0 7
8 456 2502 462 463 0 0 0 0 0 0 0 1 0 0 0 8
9 455 2593 456 462 0 0 0 0 0 0 0 0 1 0 0 9
10 456 2598 455 456 0 0 0 0 0 0 0 0 0 1 0 10
11 472 2053 456 455 0 0 0 0 0 0 0 0 0 0 1 11
12 472 2213 472 456 0 0 0 0 0 0 0 0 0 0 0 12
13 471 2238 472 472 1 0 0 0 0 0 0 0 0 0 0 13
14 465 2359 471 472 0 1 0 0 0 0 0 0 0 0 0 14
15 459 2151 465 471 0 0 1 0 0 0 0 0 0 0 0 15
16 465 2474 459 465 0 0 0 1 0 0 0 0 0 0 0 16
17 468 3079 465 459 0 0 0 0 1 0 0 0 0 0 0 17
18 467 2312 468 465 0 0 0 0 0 1 0 0 0 0 0 18
19 463 2565 467 468 0 0 0 0 0 0 1 0 0 0 0 19
20 460 1972 463 467 0 0 0 0 0 0 0 1 0 0 0 20
21 462 2484 460 463 0 0 0 0 0 0 0 0 1 0 0 21
22 461 2202 462 460 0 0 0 0 0 0 0 0 0 1 0 22
23 476 2151 461 462 0 0 0 0 0 0 0 0 0 0 1 23
24 476 1976 476 461 0 0 0 0 0 0 0 0 0 0 0 24
25 471 2012 476 476 1 0 0 0 0 0 0 0 0 0 0 25
26 453 2114 471 476 0 1 0 0 0 0 0 0 0 0 0 26
27 443 1772 453 471 0 0 1 0 0 0 0 0 0 0 0 27
28 442 1957 443 453 0 0 0 1 0 0 0 0 0 0 0 28
29 444 2070 442 443 0 0 0 0 1 0 0 0 0 0 0 29
30 438 1990 444 442 0 0 0 0 0 1 0 0 0 0 0 30
31 427 2182 438 444 0 0 0 0 0 0 1 0 0 0 0 31
32 424 2008 427 438 0 0 0 0 0 0 0 1 0 0 0 32
33 416 1916 424 427 0 0 0 0 0 0 0 0 1 0 0 33
34 406 2397 416 424 0 0 0 0 0 0 0 0 0 1 0 34
35 431 2114 406 416 0 0 0 0 0 0 0 0 0 0 1 35
36 434 1778 431 406 0 0 0 0 0 0 0 0 0 0 0 36
37 418 1641 434 431 1 0 0 0 0 0 0 0 0 0 0 37
38 412 2186 418 434 0 1 0 0 0 0 0 0 0 0 0 38
39 404 1773 412 418 0 0 1 0 0 0 0 0 0 0 0 39
40 409 1785 404 412 0 0 0 1 0 0 0 0 0 0 0 40
41 412 2217 409 404 0 0 0 0 1 0 0 0 0 0 0 41
42 406 2153 412 409 0 0 0 0 0 1 0 0 0 0 0 42
43 398 1895 406 412 0 0 0 0 0 0 1 0 0 0 0 43
44 397 2475 398 406 0 0 0 0 0 0 0 1 0 0 0 44
45 385 1793 397 398 0 0 0 0 0 0 0 0 1 0 0 45
46 390 2308 385 397 0 0 0 0 0 0 0 0 0 1 0 46
47 413 2051 390 385 0 0 0 0 0 0 0 0 0 0 1 47
48 413 1898 413 390 0 0 0 0 0 0 0 0 0 0 0 48
49 401 2142 413 413 1 0 0 0 0 0 0 0 0 0 0 49
50 397 1874 401 413 0 1 0 0 0 0 0 0 0 0 0 50
51 397 1560 397 401 0 0 1 0 0 0 0 0 0 0 0 51
52 409 1808 397 397 0 0 0 1 0 0 0 0 0 0 0 52
53 419 1575 409 397 0 0 0 0 1 0 0 0 0 0 0 53
54 424 1525 419 409 0 0 0 0 0 1 0 0 0 0 0 54
55 428 1997 424 419 0 0 0 0 0 0 1 0 0 0 0 55
56 430 1753 428 424 0 0 0 0 0 0 0 1 0 0 0 56
57 424 1623 430 428 0 0 0 0 0 0 0 0 1 0 0 57
58 433 2251 424 430 0 0 0 0 0 0 0 0 0 1 0 58
59 456 1890 433 424 0 0 0 0 0 0 0 0 0 0 1 59
> k <- length(x[1,])
> df <- as.data.frame(x)
> (mylm <- lm(df))
Call:
lm(formula = df)
Coefficients:
(Intercept) bvg Y1 Y2 M1 M2
14.182976 0.002686 1.293228 -0.348398 -0.976703 2.557338
M3 M4 M5 M6 M7 M8
1.367880 12.897522 6.564338 3.093642 1.461118 4.226885
M9 M10 M11 t
0.869449 6.674746 25.131803 -0.029334
> (mysum <- summary(mylm))
Call:
lm(formula = df)
Residuals:
Min 1Q Median 3Q Max
-11.9298 -2.8894 0.1647 3.0626 11.1305
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 14.182976 23.730103 0.598 0.5532
bvg 0.002686 0.003438 0.781 0.4388
Y1 1.293228 0.147281 8.781 3.83e-11 ***
Y2 -0.348398 0.147768 -2.358 0.0230 *
M1 -0.976703 4.551500 -0.215 0.8311
M2 2.557338 5.248527 0.487 0.6286
M3 1.367880 5.313725 0.257 0.7981
M4 12.897522 5.322429 2.423 0.0197 *
M5 6.564338 4.248698 1.545 0.1297
M6 3.093642 4.387333 0.705 0.4845
M7 1.461118 4.762780 0.307 0.7605
M8 4.226885 5.019576 0.842 0.4044
M9 0.869449 4.860084 0.179 0.8589
M10 6.674746 5.146016 1.297 0.2015
M11 25.131803 4.604406 5.458 2.23e-06 ***
t -0.029334 0.076274 -0.385 0.7024
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 5.566 on 43 degrees of freedom
Multiple R-squared: 0.9675, Adjusted R-squared: 0.9561
F-statistic: 85.21 on 15 and 43 DF, p-value: < 2.2e-16
> if (n > n25) {
+ kp3 <- k + 3
+ nmkm3 <- n - k - 3
+ gqarr <- array(NA, dim=c(nmkm3-kp3+1,3))
+ numgqtests <- 0
+ numsignificant1 <- 0
+ numsignificant5 <- 0
+ numsignificant10 <- 0
+ for (mypoint in kp3:nmkm3) {
+ j <- 0
+ numgqtests <- numgqtests + 1
+ for (myalt in c('greater', 'two.sided', 'less')) {
+ j <- j + 1
+ gqarr[mypoint-kp3+1,j] <- gqtest(mylm, point=mypoint, alternative=myalt)$p.value
+ }
+ if (gqarr[mypoint-kp3+1,2] < 0.01) numsignificant1 <- numsignificant1 + 1
+ if (gqarr[mypoint-kp3+1,2] < 0.05) numsignificant5 <- numsignificant5 + 1
+ if (gqarr[mypoint-kp3+1,2] < 0.10) numsignificant10 <- numsignificant10 + 1
+ }
+ gqarr
+ }
[,1] [,2] [,3]
[1,] 0.24207248 0.4841450 0.75792752
[2,] 0.23418024 0.4683605 0.76581976
[3,] 0.28948714 0.5789743 0.71051286
[4,] 0.20022085 0.4004417 0.79977915
[5,] 0.11594307 0.2318861 0.88405693
[6,] 0.06314551 0.1262910 0.93685449
[7,] 0.13060459 0.2612092 0.86939541
[8,] 0.74812741 0.5037452 0.25187259
[9,] 0.66131219 0.6773756 0.33868781
[10,] 0.67496337 0.6500733 0.32503663
[11,] 0.57414841 0.8517032 0.42585159
[12,] 0.52915840 0.9416832 0.47084160
[13,] 0.49285172 0.9857034 0.50714828
[14,] 0.39272384 0.7854477 0.60727616
[15,] 0.48756502 0.9751300 0.51243498
[16,] 0.65137404 0.6972519 0.34862596
[17,] 0.92578435 0.1484313 0.07421565
[18,] 0.94045392 0.1190922 0.05954608
[19,] 0.91675696 0.1664861 0.08324304
[20,] 0.93207141 0.1358572 0.06792859
[21,] 0.85688954 0.2862209 0.14311046
[22,] 0.91397097 0.1720581 0.08602903
> postscript(file="/var/www/html/rcomp/tmp/1zbkr1260820746.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/26k0l1260820746.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/3bx851260820746.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/4wgru1260820746.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/53a0v1260820746.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 = 59
Frequency = 1
1 2 3 4 5 6
6.78342972 11.13048081 -5.29940711 1.57299146 -2.60112500 4.43052877
7 8 9 10 11 12
1.93642214 -5.05981908 4.49346462 -1.09308738 -3.69837322 0.68969085
13 14 15 16 17 18
6.20293140 -2.33359677 0.85492989 0.15591181 -1.95657158 0.81459445
19 20 21 22 23 24
0.13522530 0.81631638 7.31376604 -2.33629369 -3.63698974 2.24743662
25 26 27 28 29 30
3.38272956 -11.92984447 1.74380325 -4.59235657 1.27585492 -3.94406122
31 32 33 34 35 36
-5.34181862 1.52429916 -2.79447586 -10.56194195 6.91567482 0.16474193
37 38 39 40 41 42
-9.63093365 1.13714133 -2.34959430 -0.62669893 -1.67800974 -6.14374979
43 44 45 46 47 48
-2.98424922 -0.02332903 -8.29841751 4.71248902 -0.67175532 -3.10186940
49 50 51 52 53 54
-6.73815704 1.99581909 5.05026827 3.49015223 4.95985140 4.84268779
55 56 57 58 59
6.25442040 2.74253257 -0.71433729 9.27883400 1.09144347
> postscript(file="/var/www/html/rcomp/tmp/6y2k11260820746.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 = 59
Frequency = 1
lag(myerror, k = 1) myerror
0 6.78342972 NA
1 11.13048081 6.78342972
2 -5.29940711 11.13048081
3 1.57299146 -5.29940711
4 -2.60112500 1.57299146
5 4.43052877 -2.60112500
6 1.93642214 4.43052877
7 -5.05981908 1.93642214
8 4.49346462 -5.05981908
9 -1.09308738 4.49346462
10 -3.69837322 -1.09308738
11 0.68969085 -3.69837322
12 6.20293140 0.68969085
13 -2.33359677 6.20293140
14 0.85492989 -2.33359677
15 0.15591181 0.85492989
16 -1.95657158 0.15591181
17 0.81459445 -1.95657158
18 0.13522530 0.81459445
19 0.81631638 0.13522530
20 7.31376604 0.81631638
21 -2.33629369 7.31376604
22 -3.63698974 -2.33629369
23 2.24743662 -3.63698974
24 3.38272956 2.24743662
25 -11.92984447 3.38272956
26 1.74380325 -11.92984447
27 -4.59235657 1.74380325
28 1.27585492 -4.59235657
29 -3.94406122 1.27585492
30 -5.34181862 -3.94406122
31 1.52429916 -5.34181862
32 -2.79447586 1.52429916
33 -10.56194195 -2.79447586
34 6.91567482 -10.56194195
35 0.16474193 6.91567482
36 -9.63093365 0.16474193
37 1.13714133 -9.63093365
38 -2.34959430 1.13714133
39 -0.62669893 -2.34959430
40 -1.67800974 -0.62669893
41 -6.14374979 -1.67800974
42 -2.98424922 -6.14374979
43 -0.02332903 -2.98424922
44 -8.29841751 -0.02332903
45 4.71248902 -8.29841751
46 -0.67175532 4.71248902
47 -3.10186940 -0.67175532
48 -6.73815704 -3.10186940
49 1.99581909 -6.73815704
50 5.05026827 1.99581909
51 3.49015223 5.05026827
52 4.95985140 3.49015223
53 4.84268779 4.95985140
54 6.25442040 4.84268779
55 2.74253257 6.25442040
56 -0.71433729 2.74253257
57 9.27883400 -0.71433729
58 1.09144347 9.27883400
59 NA 1.09144347
> dum1 <- dum[2:length(myerror),]
> dum1
lag(myerror, k = 1) myerror
[1,] 11.13048081 6.78342972
[2,] -5.29940711 11.13048081
[3,] 1.57299146 -5.29940711
[4,] -2.60112500 1.57299146
[5,] 4.43052877 -2.60112500
[6,] 1.93642214 4.43052877
[7,] -5.05981908 1.93642214
[8,] 4.49346462 -5.05981908
[9,] -1.09308738 4.49346462
[10,] -3.69837322 -1.09308738
[11,] 0.68969085 -3.69837322
[12,] 6.20293140 0.68969085
[13,] -2.33359677 6.20293140
[14,] 0.85492989 -2.33359677
[15,] 0.15591181 0.85492989
[16,] -1.95657158 0.15591181
[17,] 0.81459445 -1.95657158
[18,] 0.13522530 0.81459445
[19,] 0.81631638 0.13522530
[20,] 7.31376604 0.81631638
[21,] -2.33629369 7.31376604
[22,] -3.63698974 -2.33629369
[23,] 2.24743662 -3.63698974
[24,] 3.38272956 2.24743662
[25,] -11.92984447 3.38272956
[26,] 1.74380325 -11.92984447
[27,] -4.59235657 1.74380325
[28,] 1.27585492 -4.59235657
[29,] -3.94406122 1.27585492
[30,] -5.34181862 -3.94406122
[31,] 1.52429916 -5.34181862
[32,] -2.79447586 1.52429916
[33,] -10.56194195 -2.79447586
[34,] 6.91567482 -10.56194195
[35,] 0.16474193 6.91567482
[36,] -9.63093365 0.16474193
[37,] 1.13714133 -9.63093365
[38,] -2.34959430 1.13714133
[39,] -0.62669893 -2.34959430
[40,] -1.67800974 -0.62669893
[41,] -6.14374979 -1.67800974
[42,] -2.98424922 -6.14374979
[43,] -0.02332903 -2.98424922
[44,] -8.29841751 -0.02332903
[45,] 4.71248902 -8.29841751
[46,] -0.67175532 4.71248902
[47,] -3.10186940 -0.67175532
[48,] -6.73815704 -3.10186940
[49,] 1.99581909 -6.73815704
[50,] 5.05026827 1.99581909
[51,] 3.49015223 5.05026827
[52,] 4.95985140 3.49015223
[53,] 4.84268779 4.95985140
[54,] 6.25442040 4.84268779
[55,] 2.74253257 6.25442040
[56,] -0.71433729 2.74253257
[57,] 9.27883400 -0.71433729
[58,] 1.09144347 9.27883400
> z <- as.data.frame(dum1)
> z
lag(myerror, k = 1) myerror
1 11.13048081 6.78342972
2 -5.29940711 11.13048081
3 1.57299146 -5.29940711
4 -2.60112500 1.57299146
5 4.43052877 -2.60112500
6 1.93642214 4.43052877
7 -5.05981908 1.93642214
8 4.49346462 -5.05981908
9 -1.09308738 4.49346462
10 -3.69837322 -1.09308738
11 0.68969085 -3.69837322
12 6.20293140 0.68969085
13 -2.33359677 6.20293140
14 0.85492989 -2.33359677
15 0.15591181 0.85492989
16 -1.95657158 0.15591181
17 0.81459445 -1.95657158
18 0.13522530 0.81459445
19 0.81631638 0.13522530
20 7.31376604 0.81631638
21 -2.33629369 7.31376604
22 -3.63698974 -2.33629369
23 2.24743662 -3.63698974
24 3.38272956 2.24743662
25 -11.92984447 3.38272956
26 1.74380325 -11.92984447
27 -4.59235657 1.74380325
28 1.27585492 -4.59235657
29 -3.94406122 1.27585492
30 -5.34181862 -3.94406122
31 1.52429916 -5.34181862
32 -2.79447586 1.52429916
33 -10.56194195 -2.79447586
34 6.91567482 -10.56194195
35 0.16474193 6.91567482
36 -9.63093365 0.16474193
37 1.13714133 -9.63093365
38 -2.34959430 1.13714133
39 -0.62669893 -2.34959430
40 -1.67800974 -0.62669893
41 -6.14374979 -1.67800974
42 -2.98424922 -6.14374979
43 -0.02332903 -2.98424922
44 -8.29841751 -0.02332903
45 4.71248902 -8.29841751
46 -0.67175532 4.71248902
47 -3.10186940 -0.67175532
48 -6.73815704 -3.10186940
49 1.99581909 -6.73815704
50 5.05026827 1.99581909
51 3.49015223 5.05026827
52 4.95985140 3.49015223
53 4.84268779 4.95985140
54 6.25442040 4.84268779
55 2.74253257 6.25442040
56 -0.71433729 2.74253257
57 9.27883400 -0.71433729
58 1.09144347 9.27883400
> 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/75ley1260820746.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/8kgqg1260820746.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/9vwrz1260820746.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/10bvbv1260820746.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/11x5jm1260820746.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/1283p01260820746.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/13j4q21260820746.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/14h9pl1260820746.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/15u5cb1260820746.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/16cpmo1260820746.tab")
+ }
>
> try(system("convert tmp/1zbkr1260820746.ps tmp/1zbkr1260820746.png",intern=TRUE))
character(0)
> try(system("convert tmp/26k0l1260820746.ps tmp/26k0l1260820746.png",intern=TRUE))
character(0)
> try(system("convert tmp/3bx851260820746.ps tmp/3bx851260820746.png",intern=TRUE))
character(0)
> try(system("convert tmp/4wgru1260820746.ps tmp/4wgru1260820746.png",intern=TRUE))
character(0)
> try(system("convert tmp/53a0v1260820746.ps tmp/53a0v1260820746.png",intern=TRUE))
character(0)
> try(system("convert tmp/6y2k11260820746.ps tmp/6y2k11260820746.png",intern=TRUE))
character(0)
> try(system("convert tmp/75ley1260820746.ps tmp/75ley1260820746.png",intern=TRUE))
character(0)
> try(system("convert tmp/8kgqg1260820746.ps tmp/8kgqg1260820746.png",intern=TRUE))
character(0)
> try(system("convert tmp/9vwrz1260820746.ps tmp/9vwrz1260820746.png",intern=TRUE))
character(0)
> try(system("convert tmp/10bvbv1260820746.ps tmp/10bvbv1260820746.png",intern=TRUE))
character(0)
>
>
> proc.time()
user system elapsed
2.386 1.566 2.836