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(467,98.8,460,100.5,448,110.4,443,96.4,436,101.9,431,106.2,484,81,510,94.7,513,101,503,109.4,471,102.3,471,90.7,476,96.2,475,96.1,470,106,461,103.1,455,102,456,104.7,517,86,525,92.1,523,106.9,519,112.6,509,101.7,512,92,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),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 = 'No 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
1 467 98.8
2 460 100.5
3 448 110.4
4 443 96.4
5 436 101.9
6 431 106.2
7 484 81.0
8 510 94.7
9 513 101.0
10 503 109.4
11 471 102.3
12 471 90.7
13 476 96.2
14 475 96.1
15 470 106.0
16 461 103.1
17 455 102.0
18 456 104.7
19 517 86.0
20 525 92.1
21 523 106.9
22 519 112.6
23 509 101.7
24 512 92.0
25 519 97.4
26 517 97.0
27 510 105.4
28 509 102.7
29 501 98.1
30 507 104.5
31 569 87.4
32 580 89.9
33 578 109.8
34 565 111.7
35 547 98.6
36 555 96.9
37 562 95.1
38 561 97.0
39 555 112.7
40 544 102.9
41 537 97.4
42 543 111.4
43 594 87.4
44 611 96.8
45 613 114.1
46 611 110.3
47 594 103.9
48 595 101.6
49 591 94.6
50 589 95.9
51 584 104.7
52 573 102.8
53 567 98.1
54 569 113.9
55 621 80.9
56 629 95.7
57 628 113.2
58 612 105.9
59 595 108.8
60 597 102.3
61 593 99.0
62 590 100.7
63 580 115.5
64 574 100.7
65 573 109.9
66 573 114.6
67 620 85.4
68 626 100.5
69 620 114.8
70 588 116.5
71 566 112.9
72 557 102.0
> k <- length(x[1,])
> df <- as.data.frame(x)
> (mylm <- lm(df))
Call:
lm(formula = df)
Coefficients:
(Intercept) X
482.7093 0.5979
> (mysum <- summary(mylm))
Call:
lm(formula = df)
Residuals:
Min 1Q Median 3Q Max
-115.21 -36.34 13.83 47.12 89.92
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 482.7093 79.9946 6.034 6.82e-08 ***
X 0.5979 0.7851 0.762 0.449
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 55.42 on 70 degrees of freedom
Multiple R-squared: 0.008218, Adjusted R-squared: -0.00595
F-statistic: 0.58 on 1 and 70 DF, p-value: 0.4489
> 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.029156497 5.831299e-02 9.708435e-01
[2,] 0.014473262 2.894652e-02 9.855267e-01
[3,] 0.004464323 8.928646e-03 9.955357e-01
[4,] 0.028645524 5.729105e-02 9.713545e-01
[5,] 0.075174795 1.503496e-01 9.248252e-01
[6,] 0.102021186 2.040424e-01 8.979788e-01
[7,] 0.068389511 1.367790e-01 9.316105e-01
[8,] 0.045802019 9.160404e-02 9.541980e-01
[9,] 0.030506302 6.101260e-02 9.694937e-01
[10,] 0.020865936 4.173187e-02 9.791341e-01
[11,] 0.015633925 3.126785e-02 9.843661e-01
[12,] 0.013762757 2.752551e-02 9.862372e-01
[13,] 0.015720260 3.144052e-02 9.842797e-01
[14,] 0.020437323 4.087465e-02 9.795627e-01
[15,] 0.024830348 4.966070e-02 9.751697e-01
[16,] 0.040488484 8.097697e-02 9.595115e-01
[17,] 0.104572243 2.091445e-01 8.954278e-01
[18,] 0.180859054 3.617181e-01 8.191409e-01
[19,] 0.211891435 4.237829e-01 7.881086e-01
[20,] 0.235741290 4.714826e-01 7.642587e-01
[21,] 0.281508001 5.630160e-01 7.184920e-01
[22,] 0.332520659 6.650413e-01 6.674793e-01
[23,] 0.411803191 8.236064e-01 5.881968e-01
[24,] 0.509613977 9.807720e-01 4.903860e-01
[25,] 0.660719772 6.785605e-01 3.392802e-01
[26,] 0.812170338 3.756593e-01 1.878297e-01
[27,] 0.900225904 1.995482e-01 9.977410e-02
[28,] 0.949436862 1.011263e-01 5.056314e-02
[29,] 0.985071939 2.985612e-02 1.492806e-02
[30,] 0.991774609 1.645078e-02 8.225391e-03
[31,] 0.994480905 1.103819e-02 5.519095e-03
[32,] 0.996080816 7.838368e-03 3.919184e-03
[33,] 0.997046761 5.906477e-03 2.953239e-03
[34,] 0.997761509 4.476982e-03 2.238491e-03
[35,] 0.998261472 3.477055e-03 1.738528e-03
[36,] 0.999077046 1.845908e-03 9.229538e-04
[37,] 0.999789956 4.200880e-04 2.100440e-04
[38,] 0.999926728 1.465444e-04 7.327221e-05
[39,] 0.999935601 1.287981e-04 6.439905e-05
[40,] 0.999958764 8.247106e-05 4.123553e-05
[41,] 0.999983372 3.325685e-05 1.662843e-05
[42,] 0.999989188 2.162343e-05 1.081171e-05
[43,] 0.999983454 3.309118e-05 1.654559e-05
[44,] 0.999973579 5.284115e-05 2.642057e-05
[45,] 0.999957286 8.542799e-05 4.271399e-05
[46,] 0.999928497 1.430056e-04 7.150278e-05
[47,] 0.999868665 2.626691e-04 1.313346e-04
[48,] 0.999807011 3.859772e-04 1.929886e-04
[49,] 0.999824991 3.500188e-04 1.750094e-04
[50,] 0.999697336 6.053283e-04 3.026642e-04
[51,] 0.999504834 9.903326e-04 4.951663e-04
[52,] 0.999529260 9.414794e-04 4.707397e-04
[53,] 0.999802119 3.957612e-04 1.978806e-04
[54,] 0.999714387 5.712262e-04 2.856131e-04
[55,] 0.999316057 1.367887e-03 6.839435e-04
[56,] 0.998334759 3.330483e-03 1.665241e-03
[57,] 0.995922980 8.154041e-03 4.077020e-03
[58,] 0.990348435 1.930313e-02 9.651565e-03
[59,] 0.977486243 4.502751e-02 2.251376e-02
[60,] 0.960580180 7.883964e-02 3.941982e-02
[61,] 0.923352325 1.532953e-01 7.664767e-02
[62,] 0.852399084 2.952018e-01 1.476009e-01
[63,] 0.741425430 5.171491e-01 2.585746e-01
> postscript(file="/var/www/html/rcomp/tmp/1f0k61260627008.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/2xoo91260627008.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/3hllt1260627008.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/45qku1260627008.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/5hg5i1260627008.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
-74.7829765 -82.7994265 -100.7187530 -97.3479882 -107.6365029 -115.2075236
7 8 9 10 11 12
-47.1401469 -29.3315382 -30.0983823 -45.1208412 -72.8756677 -65.9398911
13 14 15 16 17 18
-64.2284058 -65.1686147 -76.0879412 -83.3539971 -88.6962941 -89.3106559
19 20 21 22 23 24
-17.1297058 -12.7769676 -23.6260618 -31.0341589 -34.5169206 -25.7171764
25 26 27 28 29 30
-21.9459000 -23.7067353 -35.7291941 -35.1148324 -40.3644382 -38.1910735
31 32 33 34 35 36
34.0332177 43.5384383 29.6399941 15.5039617 5.3366059 14.3530559
37 38 39 40 41 42
22.4292971 20.2932647 4.9060499 -0.2344147 -3.9459000 -6.3166648
43 44 45 46 47 48
59.0332177 70.4128471 62.0689735 62.3410382 49.1676735 51.5428706
49 50 51 52 53 54
51.7282530 48.9509677 38.6893441 28.8253765 25.6355618 18.1885558
55 56 57 58 59 60
89.9196442 89.0705500 77.6070940 65.9718500 47.2379058 53.1243323
61 62 63 64 65 66
51.0974412 47.0809912 28.2318970 31.0809912 24.5802029 21.7700176
67 68 69 70 71 72
86.2290413 83.2005735 68.6504352 35.6339852 15.7864676 13.3037059
> postscript(file="/var/www/html/rcomp/tmp/62dwp1260627008.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 -74.7829765 NA
1 -82.7994265 -74.7829765
2 -100.7187530 -82.7994265
3 -97.3479882 -100.7187530
4 -107.6365029 -97.3479882
5 -115.2075236 -107.6365029
6 -47.1401469 -115.2075236
7 -29.3315382 -47.1401469
8 -30.0983823 -29.3315382
9 -45.1208412 -30.0983823
10 -72.8756677 -45.1208412
11 -65.9398911 -72.8756677
12 -64.2284058 -65.9398911
13 -65.1686147 -64.2284058
14 -76.0879412 -65.1686147
15 -83.3539971 -76.0879412
16 -88.6962941 -83.3539971
17 -89.3106559 -88.6962941
18 -17.1297058 -89.3106559
19 -12.7769676 -17.1297058
20 -23.6260618 -12.7769676
21 -31.0341589 -23.6260618
22 -34.5169206 -31.0341589
23 -25.7171764 -34.5169206
24 -21.9459000 -25.7171764
25 -23.7067353 -21.9459000
26 -35.7291941 -23.7067353
27 -35.1148324 -35.7291941
28 -40.3644382 -35.1148324
29 -38.1910735 -40.3644382
30 34.0332177 -38.1910735
31 43.5384383 34.0332177
32 29.6399941 43.5384383
33 15.5039617 29.6399941
34 5.3366059 15.5039617
35 14.3530559 5.3366059
36 22.4292971 14.3530559
37 20.2932647 22.4292971
38 4.9060499 20.2932647
39 -0.2344147 4.9060499
40 -3.9459000 -0.2344147
41 -6.3166648 -3.9459000
42 59.0332177 -6.3166648
43 70.4128471 59.0332177
44 62.0689735 70.4128471
45 62.3410382 62.0689735
46 49.1676735 62.3410382
47 51.5428706 49.1676735
48 51.7282530 51.5428706
49 48.9509677 51.7282530
50 38.6893441 48.9509677
51 28.8253765 38.6893441
52 25.6355618 28.8253765
53 18.1885558 25.6355618
54 89.9196442 18.1885558
55 89.0705500 89.9196442
56 77.6070940 89.0705500
57 65.9718500 77.6070940
58 47.2379058 65.9718500
59 53.1243323 47.2379058
60 51.0974412 53.1243323
61 47.0809912 51.0974412
62 28.2318970 47.0809912
63 31.0809912 28.2318970
64 24.5802029 31.0809912
65 21.7700176 24.5802029
66 86.2290413 21.7700176
67 83.2005735 86.2290413
68 68.6504352 83.2005735
69 35.6339852 68.6504352
70 15.7864676 35.6339852
71 13.3037059 15.7864676
72 NA 13.3037059
> dum1 <- dum[2:length(myerror),]
> dum1
lag(myerror, k = 1) myerror
[1,] -82.7994265 -74.7829765
[2,] -100.7187530 -82.7994265
[3,] -97.3479882 -100.7187530
[4,] -107.6365029 -97.3479882
[5,] -115.2075236 -107.6365029
[6,] -47.1401469 -115.2075236
[7,] -29.3315382 -47.1401469
[8,] -30.0983823 -29.3315382
[9,] -45.1208412 -30.0983823
[10,] -72.8756677 -45.1208412
[11,] -65.9398911 -72.8756677
[12,] -64.2284058 -65.9398911
[13,] -65.1686147 -64.2284058
[14,] -76.0879412 -65.1686147
[15,] -83.3539971 -76.0879412
[16,] -88.6962941 -83.3539971
[17,] -89.3106559 -88.6962941
[18,] -17.1297058 -89.3106559
[19,] -12.7769676 -17.1297058
[20,] -23.6260618 -12.7769676
[21,] -31.0341589 -23.6260618
[22,] -34.5169206 -31.0341589
[23,] -25.7171764 -34.5169206
[24,] -21.9459000 -25.7171764
[25,] -23.7067353 -21.9459000
[26,] -35.7291941 -23.7067353
[27,] -35.1148324 -35.7291941
[28,] -40.3644382 -35.1148324
[29,] -38.1910735 -40.3644382
[30,] 34.0332177 -38.1910735
[31,] 43.5384383 34.0332177
[32,] 29.6399941 43.5384383
[33,] 15.5039617 29.6399941
[34,] 5.3366059 15.5039617
[35,] 14.3530559 5.3366059
[36,] 22.4292971 14.3530559
[37,] 20.2932647 22.4292971
[38,] 4.9060499 20.2932647
[39,] -0.2344147 4.9060499
[40,] -3.9459000 -0.2344147
[41,] -6.3166648 -3.9459000
[42,] 59.0332177 -6.3166648
[43,] 70.4128471 59.0332177
[44,] 62.0689735 70.4128471
[45,] 62.3410382 62.0689735
[46,] 49.1676735 62.3410382
[47,] 51.5428706 49.1676735
[48,] 51.7282530 51.5428706
[49,] 48.9509677 51.7282530
[50,] 38.6893441 48.9509677
[51,] 28.8253765 38.6893441
[52,] 25.6355618 28.8253765
[53,] 18.1885558 25.6355618
[54,] 89.9196442 18.1885558
[55,] 89.0705500 89.9196442
[56,] 77.6070940 89.0705500
[57,] 65.9718500 77.6070940
[58,] 47.2379058 65.9718500
[59,] 53.1243323 47.2379058
[60,] 51.0974412 53.1243323
[61,] 47.0809912 51.0974412
[62,] 28.2318970 47.0809912
[63,] 31.0809912 28.2318970
[64,] 24.5802029 31.0809912
[65,] 21.7700176 24.5802029
[66,] 86.2290413 21.7700176
[67,] 83.2005735 86.2290413
[68,] 68.6504352 83.2005735
[69,] 35.6339852 68.6504352
[70,] 15.7864676 35.6339852
[71,] 13.3037059 15.7864676
> z <- as.data.frame(dum1)
> z
lag(myerror, k = 1) myerror
1 -82.7994265 -74.7829765
2 -100.7187530 -82.7994265
3 -97.3479882 -100.7187530
4 -107.6365029 -97.3479882
5 -115.2075236 -107.6365029
6 -47.1401469 -115.2075236
7 -29.3315382 -47.1401469
8 -30.0983823 -29.3315382
9 -45.1208412 -30.0983823
10 -72.8756677 -45.1208412
11 -65.9398911 -72.8756677
12 -64.2284058 -65.9398911
13 -65.1686147 -64.2284058
14 -76.0879412 -65.1686147
15 -83.3539971 -76.0879412
16 -88.6962941 -83.3539971
17 -89.3106559 -88.6962941
18 -17.1297058 -89.3106559
19 -12.7769676 -17.1297058
20 -23.6260618 -12.7769676
21 -31.0341589 -23.6260618
22 -34.5169206 -31.0341589
23 -25.7171764 -34.5169206
24 -21.9459000 -25.7171764
25 -23.7067353 -21.9459000
26 -35.7291941 -23.7067353
27 -35.1148324 -35.7291941
28 -40.3644382 -35.1148324
29 -38.1910735 -40.3644382
30 34.0332177 -38.1910735
31 43.5384383 34.0332177
32 29.6399941 43.5384383
33 15.5039617 29.6399941
34 5.3366059 15.5039617
35 14.3530559 5.3366059
36 22.4292971 14.3530559
37 20.2932647 22.4292971
38 4.9060499 20.2932647
39 -0.2344147 4.9060499
40 -3.9459000 -0.2344147
41 -6.3166648 -3.9459000
42 59.0332177 -6.3166648
43 70.4128471 59.0332177
44 62.0689735 70.4128471
45 62.3410382 62.0689735
46 49.1676735 62.3410382
47 51.5428706 49.1676735
48 51.7282530 51.5428706
49 48.9509677 51.7282530
50 38.6893441 48.9509677
51 28.8253765 38.6893441
52 25.6355618 28.8253765
53 18.1885558 25.6355618
54 89.9196442 18.1885558
55 89.0705500 89.9196442
56 77.6070940 89.0705500
57 65.9718500 77.6070940
58 47.2379058 65.9718500
59 53.1243323 47.2379058
60 51.0974412 53.1243323
61 47.0809912 51.0974412
62 28.2318970 47.0809912
63 31.0809912 28.2318970
64 24.5802029 31.0809912
65 21.7700176 24.5802029
66 86.2290413 21.7700176
67 83.2005735 86.2290413
68 68.6504352 83.2005735
69 35.6339852 68.6504352
70 15.7864676 35.6339852
71 13.3037059 15.7864676
> 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/7rtg11260627008.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/84svp1260627008.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/9isge1260627008.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/10qthm1260627008.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/11y8sl1260627008.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/12rl381260627008.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/131e551260627008.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/14qatn1260627008.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/15kqt71260627008.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/16bqfk1260627008.tab")
+ }
>
> system("convert tmp/1f0k61260627008.ps tmp/1f0k61260627008.png")
> system("convert tmp/2xoo91260627008.ps tmp/2xoo91260627008.png")
> system("convert tmp/3hllt1260627008.ps tmp/3hllt1260627008.png")
> system("convert tmp/45qku1260627008.ps tmp/45qku1260627008.png")
> system("convert tmp/5hg5i1260627008.ps tmp/5hg5i1260627008.png")
> system("convert tmp/62dwp1260627008.ps tmp/62dwp1260627008.png")
> system("convert tmp/7rtg11260627008.ps tmp/7rtg11260627008.png")
> system("convert tmp/84svp1260627008.ps tmp/84svp1260627008.png")
> system("convert tmp/9isge1260627008.ps tmp/9isge1260627008.png")
> system("convert tmp/10qthm1260627008.ps tmp/10qthm1260627008.png")
>
>
> proc.time()
user system elapsed
2.573 1.566 3.218