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(519,97.4,517,97,510,105.4,509,102.7,501,98.1,507,104.5,569,87.4,580,89.9,578,109.8,565,111.7,547,98.6,555,96.9,562,95.1,561,97,555,112.7,544,102.9,537,97.4,543,111.4,594,87.4,611,96.8,613,114.1,611,110.3,594,103.9,595,101.6,591,94.6,589,95.9,584,104.7,573,102.8,567,98.1,569,113.9,621,80.9,629,95.7,628,113.2,612,105.9,595,108.8,597,102.3,593,99,590,100.7,580,115.5,574,100.7,573,109.9,573,114.6,620,85.4,626,100.5,620,114.8,588,116.5,566,112.9,557,102,561,106,549,105.3,532,118.8,526,106.1,511,109.3,499,117.2,555,92.5,565,104.2,542,112.5,527,122.4,510,113.3,514,100,517,110.7,508,112.8,493,109.8,490,117.3,469,109.1,478,115.9,528,96,534,99.8,518,116.8,506,115.7,502,99.4,516,94.3),dim=c(2,72),dimnames=list(c('Y','X'),1:72))
> y <- array(NA,dim=c(2,72),dimnames=list(c('Y','X'),1:72))
> 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 = 'Do not include Seasonal 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
Y X t
1 519 97.4 1
2 517 97.0 2
3 510 105.4 3
4 509 102.7 4
5 501 98.1 5
6 507 104.5 6
7 569 87.4 7
8 580 89.9 8
9 578 109.8 9
10 565 111.7 10
11 547 98.6 11
12 555 96.9 12
13 562 95.1 13
14 561 97.0 14
15 555 112.7 15
16 544 102.9 16
17 537 97.4 17
18 543 111.4 18
19 594 87.4 19
20 611 96.8 20
21 613 114.1 21
22 611 110.3 22
23 594 103.9 23
24 595 101.6 24
25 591 94.6 25
26 589 95.9 26
27 584 104.7 27
28 573 102.8 28
29 567 98.1 29
30 569 113.9 30
31 621 80.9 31
32 629 95.7 32
33 628 113.2 33
34 612 105.9 34
35 595 108.8 35
36 597 102.3 36
37 593 99.0 37
38 590 100.7 38
39 580 115.5 39
40 574 100.7 40
41 573 109.9 41
42 573 114.6 42
43 620 85.4 43
44 626 100.5 44
45 620 114.8 45
46 588 116.5 46
47 566 112.9 47
48 557 102.0 48
49 561 106.0 49
50 549 105.3 50
51 532 118.8 51
52 526 106.1 52
53 511 109.3 53
54 499 117.2 54
55 555 92.5 55
56 565 104.2 56
57 542 112.5 57
58 527 122.4 58
59 510 113.3 59
60 514 100.0 60
61 517 110.7 61
62 508 112.8 62
63 493 109.8 63
64 490 117.3 64
65 469 109.1 65
66 478 115.9 66
67 528 96.0 67
68 534 99.8 68
69 518 116.8 69
70 506 115.7 70
71 502 99.4 71
72 516 94.3 72
> k <- length(x[1,])
> df <- as.data.frame(x)
> (mylm <- lm(df))
Call:
lm(formula = df)
Coefficients:
(Intercept) X t
663.7472 -0.8500 -0.5146
> (mysum <- summary(mylm))
Call:
lm(formula = df)
Residuals:
Min 1Q Median 3Q Max
-76.7937 -27.1357 -0.6713 28.6920 77.4488
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 663.7472 55.0601 12.055 <2e-16 ***
X -0.8500 0.5492 -1.548 0.1263
t -0.5146 0.2343 -2.197 0.0314 *
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 38.5 on 69 degrees of freedom
Multiple R-squared: 0.1392, Adjusted R-squared: 0.1142
F-statistic: 5.577 on 2 and 69 DF, p-value: 0.005688
> 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.002627027 0.005254054 0.997372973
[2,] 0.055641218 0.111282437 0.944358782
[3,] 0.048224503 0.096449007 0.951775497
[4,] 0.186898151 0.373796302 0.813101849
[5,] 0.114909873 0.229819746 0.885090127
[6,] 0.166414816 0.332829631 0.833585184
[7,] 0.164213042 0.328426084 0.835786958
[8,] 0.141735640 0.283471280 0.858264360
[9,] 0.125948429 0.251896857 0.874051571
[10,] 0.101705975 0.203411950 0.898294025
[11,] 0.152173882 0.304347764 0.847826118
[12,] 0.318258062 0.636516124 0.681741938
[13,] 0.392569510 0.785139019 0.607430490
[14,] 0.375310676 0.750621351 0.624689324
[15,] 0.399919282 0.799838565 0.600080718
[16,] 0.465162483 0.930324965 0.534837517
[17,] 0.427774085 0.855548170 0.572225915
[18,] 0.371867061 0.743734123 0.628132939
[19,] 0.324599161 0.649198322 0.675400839
[20,] 0.320838148 0.641676297 0.679161852
[21,] 0.327675545 0.655351091 0.672324455
[22,] 0.335168413 0.670336826 0.664831587
[23,] 0.440191797 0.880383594 0.559808203
[24,] 0.654257448 0.691485104 0.345742552
[25,] 0.731182964 0.537634072 0.268817036
[26,] 0.700819947 0.598360106 0.299180053
[27,] 0.657241435 0.685517130 0.342758565
[28,] 0.662764512 0.674470976 0.337235488
[29,] 0.599520294 0.800959412 0.400479706
[30,] 0.560133732 0.879732536 0.439866268
[31,] 0.527593618 0.944812764 0.472406382
[32,] 0.521645935 0.956708130 0.478354065
[33,] 0.518299889 0.963400222 0.481700111
[34,] 0.504525582 0.990948837 0.495474418
[35,] 0.586990175 0.826019649 0.413009825
[36,] 0.600289286 0.799421428 0.399710714
[37,] 0.580658594 0.838682812 0.419341406
[38,] 0.514404094 0.971191813 0.485595906
[39,] 0.570050443 0.859899114 0.429949557
[40,] 0.805864679 0.388270641 0.194135321
[41,] 0.890353387 0.219293226 0.109646613
[42,] 0.919763375 0.160473250 0.080236625
[43,] 0.936746816 0.126506367 0.063253184
[44,] 0.946309705 0.107380591 0.053690295
[45,] 0.952363542 0.095272915 0.047636458
[46,] 0.957947988 0.084104024 0.042052012
[47,] 0.965600986 0.068798029 0.034399014
[48,] 0.977913327 0.044173346 0.022086673
[49,] 0.986351083 0.027297834 0.013648917
[50,] 0.979225194 0.041549612 0.020774806
[51,] 0.984153108 0.031693783 0.015846892
[52,] 0.986392343 0.027215314 0.013607657
[53,] 0.990617018 0.018765965 0.009382982
[54,] 0.986116848 0.027766305 0.013883152
[55,] 0.977233155 0.045533690 0.022766845
[56,] 0.972010255 0.055979490 0.027989745
[57,] 0.963768340 0.072463320 0.036231660
[58,] 0.934812501 0.130374999 0.065187499
[59,] 0.883580579 0.232838842 0.116419421
[60,] 0.909417448 0.181165104 0.090582552
[61,] 0.982739152 0.034521695 0.017260848
> postscript(file="/var/www/html/rcomp/tmp/1kzbj1258653582.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/26o5u1258653582.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/3zo1k1258653582.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/4up3b1258653582.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/5t0v51258653582.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 = 72
Frequency = 1
1 2 3 4 5 6
-61.4469551 -63.2723606 -62.6181543 -65.3984584 -76.7936787 -64.8393842
7 8 9 10 11 12
-16.8590530 -3.2195865 12.2091125 1.3386055 -27.2812398 -20.2115880
13 14 15 16 17 18
-14.2269318 -13.0974388 -5.2385546 -24.0535455 -35.2137260 -16.7997668
19 20 21 22 23 24
14.3158687 39.8200309 57.0388446 52.3235890 30.3984481 29.9581264
25 26 27 28 29 30
20.5230120 20.1425315 23.1367201 11.0363807 1.5561649 17.5000447
31 32 33 34 35 36
41.9660772 63.0600012 77.4488060 55.7587048 41.7381537 38.7280173
37 38 39 40 41 42
32.4377397 31.3972415 34.4911654 16.4263951 23.7605661 28.2699356
43 44 45 46 47 48
50.9658005 70.3147112 76.9836572 46.9431590 22.3978946 4.6479523
49 50 51 52 53 54
12.5623526 0.4819603 -4.5290584 -20.8089213 -32.5744857 -37.3452573
55 56 57 58 59 60
-1.8245910 18.6344697 3.2036804 -2.8671795 -27.0872012 -33.8770377
61 62 63 64 65 66
-21.2679329 -27.9684487 -45.0037396 -41.1144936 -68.5695551 -53.2752782
67 68 69 70 71 72
-19.6748236 -9.9304144 -10.9665875 -23.3869622 -40.7266663 -30.5468645
> postscript(file="/var/www/html/rcomp/tmp/68tlt1258653582.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 = 72
Frequency = 1
lag(myerror, k = 1) myerror
0 -61.4469551 NA
1 -63.2723606 -61.4469551
2 -62.6181543 -63.2723606
3 -65.3984584 -62.6181543
4 -76.7936787 -65.3984584
5 -64.8393842 -76.7936787
6 -16.8590530 -64.8393842
7 -3.2195865 -16.8590530
8 12.2091125 -3.2195865
9 1.3386055 12.2091125
10 -27.2812398 1.3386055
11 -20.2115880 -27.2812398
12 -14.2269318 -20.2115880
13 -13.0974388 -14.2269318
14 -5.2385546 -13.0974388
15 -24.0535455 -5.2385546
16 -35.2137260 -24.0535455
17 -16.7997668 -35.2137260
18 14.3158687 -16.7997668
19 39.8200309 14.3158687
20 57.0388446 39.8200309
21 52.3235890 57.0388446
22 30.3984481 52.3235890
23 29.9581264 30.3984481
24 20.5230120 29.9581264
25 20.1425315 20.5230120
26 23.1367201 20.1425315
27 11.0363807 23.1367201
28 1.5561649 11.0363807
29 17.5000447 1.5561649
30 41.9660772 17.5000447
31 63.0600012 41.9660772
32 77.4488060 63.0600012
33 55.7587048 77.4488060
34 41.7381537 55.7587048
35 38.7280173 41.7381537
36 32.4377397 38.7280173
37 31.3972415 32.4377397
38 34.4911654 31.3972415
39 16.4263951 34.4911654
40 23.7605661 16.4263951
41 28.2699356 23.7605661
42 50.9658005 28.2699356
43 70.3147112 50.9658005
44 76.9836572 70.3147112
45 46.9431590 76.9836572
46 22.3978946 46.9431590
47 4.6479523 22.3978946
48 12.5623526 4.6479523
49 0.4819603 12.5623526
50 -4.5290584 0.4819603
51 -20.8089213 -4.5290584
52 -32.5744857 -20.8089213
53 -37.3452573 -32.5744857
54 -1.8245910 -37.3452573
55 18.6344697 -1.8245910
56 3.2036804 18.6344697
57 -2.8671795 3.2036804
58 -27.0872012 -2.8671795
59 -33.8770377 -27.0872012
60 -21.2679329 -33.8770377
61 -27.9684487 -21.2679329
62 -45.0037396 -27.9684487
63 -41.1144936 -45.0037396
64 -68.5695551 -41.1144936
65 -53.2752782 -68.5695551
66 -19.6748236 -53.2752782
67 -9.9304144 -19.6748236
68 -10.9665875 -9.9304144
69 -23.3869622 -10.9665875
70 -40.7266663 -23.3869622
71 -30.5468645 -40.7266663
72 NA -30.5468645
> dum1 <- dum[2:length(myerror),]
> dum1
lag(myerror, k = 1) myerror
[1,] -63.2723606 -61.4469551
[2,] -62.6181543 -63.2723606
[3,] -65.3984584 -62.6181543
[4,] -76.7936787 -65.3984584
[5,] -64.8393842 -76.7936787
[6,] -16.8590530 -64.8393842
[7,] -3.2195865 -16.8590530
[8,] 12.2091125 -3.2195865
[9,] 1.3386055 12.2091125
[10,] -27.2812398 1.3386055
[11,] -20.2115880 -27.2812398
[12,] -14.2269318 -20.2115880
[13,] -13.0974388 -14.2269318
[14,] -5.2385546 -13.0974388
[15,] -24.0535455 -5.2385546
[16,] -35.2137260 -24.0535455
[17,] -16.7997668 -35.2137260
[18,] 14.3158687 -16.7997668
[19,] 39.8200309 14.3158687
[20,] 57.0388446 39.8200309
[21,] 52.3235890 57.0388446
[22,] 30.3984481 52.3235890
[23,] 29.9581264 30.3984481
[24,] 20.5230120 29.9581264
[25,] 20.1425315 20.5230120
[26,] 23.1367201 20.1425315
[27,] 11.0363807 23.1367201
[28,] 1.5561649 11.0363807
[29,] 17.5000447 1.5561649
[30,] 41.9660772 17.5000447
[31,] 63.0600012 41.9660772
[32,] 77.4488060 63.0600012
[33,] 55.7587048 77.4488060
[34,] 41.7381537 55.7587048
[35,] 38.7280173 41.7381537
[36,] 32.4377397 38.7280173
[37,] 31.3972415 32.4377397
[38,] 34.4911654 31.3972415
[39,] 16.4263951 34.4911654
[40,] 23.7605661 16.4263951
[41,] 28.2699356 23.7605661
[42,] 50.9658005 28.2699356
[43,] 70.3147112 50.9658005
[44,] 76.9836572 70.3147112
[45,] 46.9431590 76.9836572
[46,] 22.3978946 46.9431590
[47,] 4.6479523 22.3978946
[48,] 12.5623526 4.6479523
[49,] 0.4819603 12.5623526
[50,] -4.5290584 0.4819603
[51,] -20.8089213 -4.5290584
[52,] -32.5744857 -20.8089213
[53,] -37.3452573 -32.5744857
[54,] -1.8245910 -37.3452573
[55,] 18.6344697 -1.8245910
[56,] 3.2036804 18.6344697
[57,] -2.8671795 3.2036804
[58,] -27.0872012 -2.8671795
[59,] -33.8770377 -27.0872012
[60,] -21.2679329 -33.8770377
[61,] -27.9684487 -21.2679329
[62,] -45.0037396 -27.9684487
[63,] -41.1144936 -45.0037396
[64,] -68.5695551 -41.1144936
[65,] -53.2752782 -68.5695551
[66,] -19.6748236 -53.2752782
[67,] -9.9304144 -19.6748236
[68,] -10.9665875 -9.9304144
[69,] -23.3869622 -10.9665875
[70,] -40.7266663 -23.3869622
[71,] -30.5468645 -40.7266663
> z <- as.data.frame(dum1)
> z
lag(myerror, k = 1) myerror
1 -63.2723606 -61.4469551
2 -62.6181543 -63.2723606
3 -65.3984584 -62.6181543
4 -76.7936787 -65.3984584
5 -64.8393842 -76.7936787
6 -16.8590530 -64.8393842
7 -3.2195865 -16.8590530
8 12.2091125 -3.2195865
9 1.3386055 12.2091125
10 -27.2812398 1.3386055
11 -20.2115880 -27.2812398
12 -14.2269318 -20.2115880
13 -13.0974388 -14.2269318
14 -5.2385546 -13.0974388
15 -24.0535455 -5.2385546
16 -35.2137260 -24.0535455
17 -16.7997668 -35.2137260
18 14.3158687 -16.7997668
19 39.8200309 14.3158687
20 57.0388446 39.8200309
21 52.3235890 57.0388446
22 30.3984481 52.3235890
23 29.9581264 30.3984481
24 20.5230120 29.9581264
25 20.1425315 20.5230120
26 23.1367201 20.1425315
27 11.0363807 23.1367201
28 1.5561649 11.0363807
29 17.5000447 1.5561649
30 41.9660772 17.5000447
31 63.0600012 41.9660772
32 77.4488060 63.0600012
33 55.7587048 77.4488060
34 41.7381537 55.7587048
35 38.7280173 41.7381537
36 32.4377397 38.7280173
37 31.3972415 32.4377397
38 34.4911654 31.3972415
39 16.4263951 34.4911654
40 23.7605661 16.4263951
41 28.2699356 23.7605661
42 50.9658005 28.2699356
43 70.3147112 50.9658005
44 76.9836572 70.3147112
45 46.9431590 76.9836572
46 22.3978946 46.9431590
47 4.6479523 22.3978946
48 12.5623526 4.6479523
49 0.4819603 12.5623526
50 -4.5290584 0.4819603
51 -20.8089213 -4.5290584
52 -32.5744857 -20.8089213
53 -37.3452573 -32.5744857
54 -1.8245910 -37.3452573
55 18.6344697 -1.8245910
56 3.2036804 18.6344697
57 -2.8671795 3.2036804
58 -27.0872012 -2.8671795
59 -33.8770377 -27.0872012
60 -21.2679329 -33.8770377
61 -27.9684487 -21.2679329
62 -45.0037396 -27.9684487
63 -41.1144936 -45.0037396
64 -68.5695551 -41.1144936
65 -53.2752782 -68.5695551
66 -19.6748236 -53.2752782
67 -9.9304144 -19.6748236
68 -10.9665875 -9.9304144
69 -23.3869622 -10.9665875
70 -40.7266663 -23.3869622
71 -30.5468645 -40.7266663
> 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/7k7tv1258653582.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/84xxz1258653582.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/9q8p31258653582.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/104jcz1258653582.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/11g4p41258653582.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/12g7vq1258653582.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/13loxx1258653583.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/14vgfx1258653583.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/156l991258653583.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/166ua11258653583.tab")
+ }
>
> system("convert tmp/1kzbj1258653582.ps tmp/1kzbj1258653582.png")
> system("convert tmp/26o5u1258653582.ps tmp/26o5u1258653582.png")
> system("convert tmp/3zo1k1258653582.ps tmp/3zo1k1258653582.png")
> system("convert tmp/4up3b1258653582.ps tmp/4up3b1258653582.png")
> system("convert tmp/5t0v51258653582.ps tmp/5t0v51258653582.png")
> system("convert tmp/68tlt1258653582.ps tmp/68tlt1258653582.png")
> system("convert tmp/7k7tv1258653582.ps tmp/7k7tv1258653582.png")
> system("convert tmp/84xxz1258653582.ps tmp/84xxz1258653582.png")
> system("convert tmp/9q8p31258653582.ps tmp/9q8p31258653582.png")
> system("convert tmp/104jcz1258653582.ps tmp/104jcz1258653582.png")
>
>
> proc.time()
user system elapsed
2.592 1.559 7.407