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(2649.2,31077,2579.4,31293,2504.6,30236,2462.3,30160,2467.4,32436,2446.7,30695,2656.3,27525,2626.2,26434,2482.6,25739,2539.9,25204,2502.7,24977,2466.9,24320,2513.2,22680,2443.3,22052,2293.4,21467,2070.8,21383,2029.6,21777,2052,21928,1864.4,21814,1670.1,22937,1811,23595,1905.4,20830,1862.8,19650,2014.5,19195,2197.8,19644,2962.3,18483,3047,18079,3032.6,19178,3504.4,18391,3801.1,18441,3857.6,18584,3674.4,20108,3721,20148,3844.5,19394,4116.7,17745,4105.2,17696,4435.2,17032,4296.5,16438,4202.5,15683,4562.8,15594,4621.4,15713,4697,15937,4591.3,16171,4357,15928,4502.6,16348,4443.9,15579,4290.9,15305,4199.8,15648,4138.5,14954,3970.1,15137,3862.3,15839,3701.6,16050,3570.12,15168,3801.06,17064,3895.51,16005,3917.96,14886,3813.06,14931,3667.03,14544,3494.17,13812,3364,13031,3295.3,12574,3277.0,11964,3257.2,11451,3161.7,11346,3097.3,11353,3061.3,10702,3119.3,10646,3106.22,10556,3080.58,10463,2981.85,10407),dim=c(2,70),dimnames=list(c('Bel20','Goudprijs'),1:70))
> y <- array(NA,dim=c(2,70),dimnames=list(c('Bel20','Goudprijs'),1:70))
> 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
Bel20 Goudprijs
1 2649.20 31077
2 2579.40 31293
3 2504.60 30236
4 2462.30 30160
5 2467.40 32436
6 2446.70 30695
7 2656.30 27525
8 2626.20 26434
9 2482.60 25739
10 2539.90 25204
11 2502.70 24977
12 2466.90 24320
13 2513.20 22680
14 2443.30 22052
15 2293.40 21467
16 2070.80 21383
17 2029.60 21777
18 2052.00 21928
19 1864.40 21814
20 1670.10 22937
21 1811.00 23595
22 1905.40 20830
23 1862.80 19650
24 2014.50 19195
25 2197.80 19644
26 2962.30 18483
27 3047.00 18079
28 3032.60 19178
29 3504.40 18391
30 3801.10 18441
31 3857.60 18584
32 3674.40 20108
33 3721.00 20148
34 3844.50 19394
35 4116.70 17745
36 4105.20 17696
37 4435.20 17032
38 4296.50 16438
39 4202.50 15683
40 4562.80 15594
41 4621.40 15713
42 4697.00 15937
43 4591.30 16171
44 4357.00 15928
45 4502.60 16348
46 4443.90 15579
47 4290.90 15305
48 4199.80 15648
49 4138.50 14954
50 3970.10 15137
51 3862.30 15839
52 3701.60 16050
53 3570.12 15168
54 3801.06 17064
55 3895.51 16005
56 3917.96 14886
57 3813.06 14931
58 3667.03 14544
59 3494.17 13812
60 3364.00 13031
61 3295.30 12574
62 3277.00 11964
63 3257.20 11451
64 3161.70 11346
65 3097.30 11353
66 3061.30 10702
67 3119.30 10646
68 3106.22 10556
69 3080.58 10463
70 2981.85 10407
> k <- length(x[1,])
> df <- as.data.frame(x)
> (mylm <- lm(df))
Call:
lm(formula = df)
Coefficients:
(Intercept) Goudprijs
4765.30342 -0.08168
> (mysum <- summary(mylm))
Call:
lm(formula = df)
Residuals:
Min 1Q Median 3Q Max
-1297.6 -650.9 114.4 583.8 1233.4
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 4765.30342 301.41618 15.810 < 2e-16 ***
Goudprijs -0.08168 0.01544 -5.289 1.41e-06 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 719.6 on 68 degrees of freedom
Multiple R-squared: 0.2915, Adjusted R-squared: 0.2811
F-statistic: 27.98 on 1 and 68 DF, p-value: 1.411e-06
> 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,] 2.717823e-03 5.435645e-03 9.972822e-01
[2,] 4.150892e-04 8.301784e-04 9.995849e-01
[3,] 6.685188e-05 1.337038e-04 9.999331e-01
[4,] 6.899003e-06 1.379801e-05 9.999931e-01
[5,] 2.221268e-06 4.442535e-06 9.999978e-01
[6,] 2.537843e-07 5.075686e-07 9.999997e-01
[7,] 3.247324e-08 6.494648e-08 1.000000e+00
[8,] 4.914334e-09 9.828669e-09 1.000000e+00
[9,] 4.816551e-10 9.633102e-10 1.000000e+00
[10,] 6.984086e-11 1.396817e-10 1.000000e+00
[11,] 8.251954e-11 1.650391e-10 1.000000e+00
[12,] 1.876461e-09 3.752922e-09 1.000000e+00
[13,] 6.245995e-09 1.249199e-08 1.000000e+00
[14,] 6.358380e-09 1.271676e-08 1.000000e+00
[15,] 3.114730e-08 6.229459e-08 1.000000e+00
[16,] 7.919740e-07 1.583948e-06 9.999992e-01
[17,] 2.767529e-06 5.535057e-06 9.999972e-01
[18,] 3.897376e-06 7.794752e-06 9.999961e-01
[19,] 8.867793e-06 1.773559e-05 9.999911e-01
[20,] 2.513404e-05 5.026807e-05 9.999749e-01
[21,] 1.546901e-04 3.093801e-04 9.998453e-01
[22,] 1.087284e-02 2.174568e-02 9.891272e-01
[23,] 9.323990e-02 1.864798e-01 9.067601e-01
[24,] 3.291690e-01 6.583380e-01 6.708310e-01
[25,] 7.035013e-01 5.929974e-01 2.964987e-01
[26,] 9.144343e-01 1.711314e-01 8.556570e-02
[27,] 9.736949e-01 5.261022e-02 2.630511e-02
[28,] 9.952226e-01 9.554722e-03 4.777361e-03
[29,] 9.996883e-01 6.234070e-04 3.117035e-04
[30,] 9.999882e-01 2.365156e-05 1.182578e-05
[31,] 9.999966e-01 6.735118e-06 3.367559e-06
[32,] 9.999989e-01 2.189696e-06 1.094848e-06
[33,] 9.999993e-01 1.455357e-06 7.276784e-07
[34,] 9.999991e-01 1.749397e-06 8.746986e-07
[35,] 9.999986e-01 2.858171e-06 1.429086e-06
[36,] 9.999993e-01 1.301392e-06 6.506958e-07
[37,] 9.999998e-01 4.277327e-07 2.138663e-07
[38,] 1.000000e+00 7.463528e-08 3.731764e-08
[39,] 1.000000e+00 2.810246e-08 1.405123e-08
[40,] 1.000000e+00 3.421635e-08 1.710818e-08
[41,] 1.000000e+00 1.864421e-08 9.322104e-09
[42,] 1.000000e+00 2.502692e-09 1.251346e-09
[43,] 1.000000e+00 3.984954e-10 1.992477e-10
[44,] 1.000000e+00 1.656132e-10 8.280659e-11
[45,] 1.000000e+00 5.731075e-12 2.865538e-12
[46,] 1.000000e+00 2.419965e-12 1.209983e-12
[47,] 1.000000e+00 1.414159e-11 7.070797e-12
[48,] 1.000000e+00 4.974251e-11 2.487125e-11
[49,] 1.000000e+00 1.095878e-10 5.479392e-11
[50,] 1.000000e+00 1.246014e-11 6.230068e-12
[51,] 1.000000e+00 1.019192e-10 5.095958e-11
[52,] 1.000000e+00 3.128306e-11 1.564153e-11
[53,] 1.000000e+00 3.845280e-11 1.922640e-11
[54,] 1.000000e+00 1.876921e-10 9.384605e-11
[55,] 1.000000e+00 2.434096e-09 1.217048e-09
[56,] 1.000000e+00 3.077682e-08 1.538841e-08
[57,] 9.999999e-01 2.884051e-07 1.442026e-07
[58,] 9.999983e-01 3.391767e-06 1.695884e-06
[59,] 9.999945e-01 1.091514e-05 5.457570e-06
[60,] 9.999395e-01 1.209007e-04 6.045035e-05
[61,] 9.994575e-01 1.084970e-03 5.424848e-04
> postscript(file="/var/www/html/rcomp/tmp/1962x1292275047.ps",horizontal=F,onefile=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/2962x1292275047.ps",horizontal=F,onefile=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/3ky101292275047.ps",horizontal=F,onefile=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/4ky101292275047.ps",horizontal=F,onefile=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/5ky101292275047.ps",horizontal=F,onefile=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 = 70
Frequency = 1
1 2 3 4 5 6
422.11190 369.95371 208.82319 160.31588 351.30828 188.41203
7 8 9 10 11 12
139.10214 19.89467 -180.46948 -166.86563 -222.60587 -312.06637
13 14 15 16 17 18
-399.71344 -520.90537 -718.58527 -948.04597 -957.06601 -922.33307
19 20 21 22 23 24
-1119.24403 -1221.82296 -1027.18078 -1158.61227 -1297.58882 -1183.05097
25 26 27 28 29 30
-963.07887 -293.40360 -241.70032 -166.33944 241.18230 541.96605
31 32 33 34 35 36
610.14558 551.41835 601.28535 663.20237 800.72022 785.21814
37 38 39 40 41 42
1060.98591 873.77094 718.10628 1071.13720 1139.45653 1233.35174
43 44 45 46 47 48
1146.76370 892.61667 1072.52018 951.01208 775.63311 712.54765
49 50 51 52 53 54
594.56517 441.11171 390.64759 247.18102 43.66363 429.45952
55 56 57 58 59 60
437.41564 368.47127 267.24665 89.60841 -143.03773 -336.99593
61 62 63 64 65 66
-443.02143 -511.14321 -572.84250 -676.91838 -740.74666 -829.91711
67 68 69 70
-776.49091 -796.92167 -830.15745 -933.46125
> postscript(file="/var/www/html/rcomp/tmp/6ky101292275047.ps",horizontal=F,onefile=F,pagecentre=F,paper="special",width=8.3333333333333,height=5.5555555555556)
> dum <- cbind(lag(myerror,k=1),myerror)
> dum
Time Series:
Start = 0
End = 70
Frequency = 1
lag(myerror, k = 1) myerror
0 422.11190 NA
1 369.95371 422.11190
2 208.82319 369.95371
3 160.31588 208.82319
4 351.30828 160.31588
5 188.41203 351.30828
6 139.10214 188.41203
7 19.89467 139.10214
8 -180.46948 19.89467
9 -166.86563 -180.46948
10 -222.60587 -166.86563
11 -312.06637 -222.60587
12 -399.71344 -312.06637
13 -520.90537 -399.71344
14 -718.58527 -520.90537
15 -948.04597 -718.58527
16 -957.06601 -948.04597
17 -922.33307 -957.06601
18 -1119.24403 -922.33307
19 -1221.82296 -1119.24403
20 -1027.18078 -1221.82296
21 -1158.61227 -1027.18078
22 -1297.58882 -1158.61227
23 -1183.05097 -1297.58882
24 -963.07887 -1183.05097
25 -293.40360 -963.07887
26 -241.70032 -293.40360
27 -166.33944 -241.70032
28 241.18230 -166.33944
29 541.96605 241.18230
30 610.14558 541.96605
31 551.41835 610.14558
32 601.28535 551.41835
33 663.20237 601.28535
34 800.72022 663.20237
35 785.21814 800.72022
36 1060.98591 785.21814
37 873.77094 1060.98591
38 718.10628 873.77094
39 1071.13720 718.10628
40 1139.45653 1071.13720
41 1233.35174 1139.45653
42 1146.76370 1233.35174
43 892.61667 1146.76370
44 1072.52018 892.61667
45 951.01208 1072.52018
46 775.63311 951.01208
47 712.54765 775.63311
48 594.56517 712.54765
49 441.11171 594.56517
50 390.64759 441.11171
51 247.18102 390.64759
52 43.66363 247.18102
53 429.45952 43.66363
54 437.41564 429.45952
55 368.47127 437.41564
56 267.24665 368.47127
57 89.60841 267.24665
58 -143.03773 89.60841
59 -336.99593 -143.03773
60 -443.02143 -336.99593
61 -511.14321 -443.02143
62 -572.84250 -511.14321
63 -676.91838 -572.84250
64 -740.74666 -676.91838
65 -829.91711 -740.74666
66 -776.49091 -829.91711
67 -796.92167 -776.49091
68 -830.15745 -796.92167
69 -933.46125 -830.15745
70 NA -933.46125
> dum1 <- dum[2:length(myerror),]
> dum1
lag(myerror, k = 1) myerror
[1,] 369.95371 422.11190
[2,] 208.82319 369.95371
[3,] 160.31588 208.82319
[4,] 351.30828 160.31588
[5,] 188.41203 351.30828
[6,] 139.10214 188.41203
[7,] 19.89467 139.10214
[8,] -180.46948 19.89467
[9,] -166.86563 -180.46948
[10,] -222.60587 -166.86563
[11,] -312.06637 -222.60587
[12,] -399.71344 -312.06637
[13,] -520.90537 -399.71344
[14,] -718.58527 -520.90537
[15,] -948.04597 -718.58527
[16,] -957.06601 -948.04597
[17,] -922.33307 -957.06601
[18,] -1119.24403 -922.33307
[19,] -1221.82296 -1119.24403
[20,] -1027.18078 -1221.82296
[21,] -1158.61227 -1027.18078
[22,] -1297.58882 -1158.61227
[23,] -1183.05097 -1297.58882
[24,] -963.07887 -1183.05097
[25,] -293.40360 -963.07887
[26,] -241.70032 -293.40360
[27,] -166.33944 -241.70032
[28,] 241.18230 -166.33944
[29,] 541.96605 241.18230
[30,] 610.14558 541.96605
[31,] 551.41835 610.14558
[32,] 601.28535 551.41835
[33,] 663.20237 601.28535
[34,] 800.72022 663.20237
[35,] 785.21814 800.72022
[36,] 1060.98591 785.21814
[37,] 873.77094 1060.98591
[38,] 718.10628 873.77094
[39,] 1071.13720 718.10628
[40,] 1139.45653 1071.13720
[41,] 1233.35174 1139.45653
[42,] 1146.76370 1233.35174
[43,] 892.61667 1146.76370
[44,] 1072.52018 892.61667
[45,] 951.01208 1072.52018
[46,] 775.63311 951.01208
[47,] 712.54765 775.63311
[48,] 594.56517 712.54765
[49,] 441.11171 594.56517
[50,] 390.64759 441.11171
[51,] 247.18102 390.64759
[52,] 43.66363 247.18102
[53,] 429.45952 43.66363
[54,] 437.41564 429.45952
[55,] 368.47127 437.41564
[56,] 267.24665 368.47127
[57,] 89.60841 267.24665
[58,] -143.03773 89.60841
[59,] -336.99593 -143.03773
[60,] -443.02143 -336.99593
[61,] -511.14321 -443.02143
[62,] -572.84250 -511.14321
[63,] -676.91838 -572.84250
[64,] -740.74666 -676.91838
[65,] -829.91711 -740.74666
[66,] -776.49091 -829.91711
[67,] -796.92167 -776.49091
[68,] -830.15745 -796.92167
[69,] -933.46125 -830.15745
> z <- as.data.frame(dum1)
> z
lag(myerror, k = 1) myerror
1 369.95371 422.11190
2 208.82319 369.95371
3 160.31588 208.82319
4 351.30828 160.31588
5 188.41203 351.30828
6 139.10214 188.41203
7 19.89467 139.10214
8 -180.46948 19.89467
9 -166.86563 -180.46948
10 -222.60587 -166.86563
11 -312.06637 -222.60587
12 -399.71344 -312.06637
13 -520.90537 -399.71344
14 -718.58527 -520.90537
15 -948.04597 -718.58527
16 -957.06601 -948.04597
17 -922.33307 -957.06601
18 -1119.24403 -922.33307
19 -1221.82296 -1119.24403
20 -1027.18078 -1221.82296
21 -1158.61227 -1027.18078
22 -1297.58882 -1158.61227
23 -1183.05097 -1297.58882
24 -963.07887 -1183.05097
25 -293.40360 -963.07887
26 -241.70032 -293.40360
27 -166.33944 -241.70032
28 241.18230 -166.33944
29 541.96605 241.18230
30 610.14558 541.96605
31 551.41835 610.14558
32 601.28535 551.41835
33 663.20237 601.28535
34 800.72022 663.20237
35 785.21814 800.72022
36 1060.98591 785.21814
37 873.77094 1060.98591
38 718.10628 873.77094
39 1071.13720 718.10628
40 1139.45653 1071.13720
41 1233.35174 1139.45653
42 1146.76370 1233.35174
43 892.61667 1146.76370
44 1072.52018 892.61667
45 951.01208 1072.52018
46 775.63311 951.01208
47 712.54765 775.63311
48 594.56517 712.54765
49 441.11171 594.56517
50 390.64759 441.11171
51 247.18102 390.64759
52 43.66363 247.18102
53 429.45952 43.66363
54 437.41564 429.45952
55 368.47127 437.41564
56 267.24665 368.47127
57 89.60841 267.24665
58 -143.03773 89.60841
59 -336.99593 -143.03773
60 -443.02143 -336.99593
61 -511.14321 -443.02143
62 -572.84250 -511.14321
63 -676.91838 -572.84250
64 -740.74666 -676.91838
65 -829.91711 -740.74666
66 -776.49091 -829.91711
67 -796.92167 -776.49091
68 -830.15745 -796.92167
69 -933.46125 -830.15745
> 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/7c7031292275047.ps",horizontal=F,onefile=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/85yzo1292275047.ps",horizontal=F,onefile=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/95yzo1292275047.ps",horizontal=F,onefile=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/105yzo1292275047.ps",horizontal=F,onefile=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/11jqxf1292275047.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/12cze01292275047.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/131icc1292275047.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/14t9bw1292275047.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/15xar21292275047.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/16bk7b1292275047.tab")
+ }
>
> try(system("convert tmp/1962x1292275047.ps tmp/1962x1292275047.png",intern=TRUE))
character(0)
> try(system("convert tmp/2962x1292275047.ps tmp/2962x1292275047.png",intern=TRUE))
character(0)
> try(system("convert tmp/3ky101292275047.ps tmp/3ky101292275047.png",intern=TRUE))
character(0)
> try(system("convert tmp/4ky101292275047.ps tmp/4ky101292275047.png",intern=TRUE))
character(0)
> try(system("convert tmp/5ky101292275047.ps tmp/5ky101292275047.png",intern=TRUE))
character(0)
> try(system("convert tmp/6ky101292275047.ps tmp/6ky101292275047.png",intern=TRUE))
character(0)
> try(system("convert tmp/7c7031292275047.ps tmp/7c7031292275047.png",intern=TRUE))
character(0)
> try(system("convert tmp/85yzo1292275047.ps tmp/85yzo1292275047.png",intern=TRUE))
character(0)
> try(system("convert tmp/95yzo1292275047.ps tmp/95yzo1292275047.png",intern=TRUE))
character(0)
> try(system("convert tmp/105yzo1292275047.ps tmp/105yzo1292275047.png",intern=TRUE))
character(0)
>
>
> proc.time()
user system elapsed
2.572 1.628 7.425