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(416.25,1111.92,398.35,1131.13,400.00,1144.94,427.25,1113.89,391.25,1107.30,397.20,1120.68,394.80,1140.84,391.50,1101.72,407.65,1104.24,418.10,1114.58,429.10,1130.20,452.85,1173.78,427.75,1211.92,420.90,1181.27,433.45,1203.60,427.15,1180.59,427.90,1156.85,415.35,1191.50,432.60,1191.33,431.65,1234.18,439.60,1220.33,466.10,1228.81,459.50,1207.01,499.75,1249.48,530.00,1248.29,568.25,1280.08,564.25,1280.66,587.00,1302.88,661.00,1310.61,625.00,1270.05,622.95,1270.06,637.25,1278.53,621.05,1303.80,600.60,1335.83,614.10,1377.76,648.75,1400.63,639.75,1418.03,660.20,1437.90,670.40,1406.80,658.25,1420.83,673.60,1482.37,666.50,1530.63,654.75,1504.66,665.75,1455.18,672.00,1473.96,742.50,1527.29,790.25,1545.79,784.25,1479.63,846.75,1467.97,914.75,1378.60,988.50,1330.45,887.75,1326.41,853.00,1385.97,888.25,1399.62,937.50,1276.69,912.50,1269.42,822.25,1287.83,880.00,1164.17,729.50,968.67,778.00,888.61),dim=c(2,60),dimnames=list(c('Gold','S&P500'),1:60))
> y <- array(NA,dim=c(2,60),dimnames=list(c('Gold','S&P500'),1:60))
> 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 = 'Include Monthly Dummies'
> par1 = '2'
> #'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
S&P500 Gold M1 M2 M3 M4 M5 M6 M7 M8 M9 M10 M11
1 1111.92 416.25 1 0 0 0 0 0 0 0 0 0 0
2 1131.13 398.35 0 1 0 0 0 0 0 0 0 0 0
3 1144.94 400.00 0 0 1 0 0 0 0 0 0 0 0
4 1113.89 427.25 0 0 0 1 0 0 0 0 0 0 0
5 1107.30 391.25 0 0 0 0 1 0 0 0 0 0 0
6 1120.68 397.20 0 0 0 0 0 1 0 0 0 0 0
7 1140.84 394.80 0 0 0 0 0 0 1 0 0 0 0
8 1101.72 391.50 0 0 0 0 0 0 0 1 0 0 0
9 1104.24 407.65 0 0 0 0 0 0 0 0 1 0 0
10 1114.58 418.10 0 0 0 0 0 0 0 0 0 1 0
11 1130.20 429.10 0 0 0 0 0 0 0 0 0 0 1
12 1173.78 452.85 0 0 0 0 0 0 0 0 0 0 0
13 1211.92 427.75 1 0 0 0 0 0 0 0 0 0 0
14 1181.27 420.90 0 1 0 0 0 0 0 0 0 0 0
15 1203.60 433.45 0 0 1 0 0 0 0 0 0 0 0
16 1180.59 427.15 0 0 0 1 0 0 0 0 0 0 0
17 1156.85 427.90 0 0 0 0 1 0 0 0 0 0 0
18 1191.50 415.35 0 0 0 0 0 1 0 0 0 0 0
19 1191.33 432.60 0 0 0 0 0 0 1 0 0 0 0
20 1234.18 431.65 0 0 0 0 0 0 0 1 0 0 0
21 1220.33 439.60 0 0 0 0 0 0 0 0 1 0 0
22 1228.81 466.10 0 0 0 0 0 0 0 0 0 1 0
23 1207.01 459.50 0 0 0 0 0 0 0 0 0 0 1
24 1249.48 499.75 0 0 0 0 0 0 0 0 0 0 0
25 1248.29 530.00 1 0 0 0 0 0 0 0 0 0 0
26 1280.08 568.25 0 1 0 0 0 0 0 0 0 0 0
27 1280.66 564.25 0 0 1 0 0 0 0 0 0 0 0
28 1302.88 587.00 0 0 0 1 0 0 0 0 0 0 0
29 1310.61 661.00 0 0 0 0 1 0 0 0 0 0 0
30 1270.05 625.00 0 0 0 0 0 1 0 0 0 0 0
31 1270.06 622.95 0 0 0 0 0 0 1 0 0 0 0
32 1278.53 637.25 0 0 0 0 0 0 0 1 0 0 0
33 1303.80 621.05 0 0 0 0 0 0 0 0 1 0 0
34 1335.83 600.60 0 0 0 0 0 0 0 0 0 1 0
35 1377.76 614.10 0 0 0 0 0 0 0 0 0 0 1
36 1400.63 648.75 0 0 0 0 0 0 0 0 0 0 0
37 1418.03 639.75 1 0 0 0 0 0 0 0 0 0 0
38 1437.90 660.20 0 1 0 0 0 0 0 0 0 0 0
39 1406.80 670.40 0 0 1 0 0 0 0 0 0 0 0
40 1420.83 658.25 0 0 0 1 0 0 0 0 0 0 0
41 1482.37 673.60 0 0 0 0 1 0 0 0 0 0 0
42 1530.63 666.50 0 0 0 0 0 1 0 0 0 0 0
43 1504.66 654.75 0 0 0 0 0 0 1 0 0 0 0
44 1455.18 665.75 0 0 0 0 0 0 0 1 0 0 0
45 1473.96 672.00 0 0 0 0 0 0 0 0 1 0 0
46 1527.29 742.50 0 0 0 0 0 0 0 0 0 1 0
47 1545.79 790.25 0 0 0 0 0 0 0 0 0 0 1
48 1479.63 784.25 0 0 0 0 0 0 0 0 0 0 0
49 1467.97 846.75 1 0 0 0 0 0 0 0 0 0 0
50 1378.60 914.75 0 1 0 0 0 0 0 0 0 0 0
51 1330.45 988.50 0 0 1 0 0 0 0 0 0 0 0
52 1326.41 887.75 0 0 0 1 0 0 0 0 0 0 0
53 1385.97 853.00 0 0 0 0 1 0 0 0 0 0 0
54 1399.62 888.25 0 0 0 0 0 1 0 0 0 0 0
55 1276.69 937.50 0 0 0 0 0 0 1 0 0 0 0
56 1269.42 912.50 0 0 0 0 0 0 0 1 0 0 0
57 1287.83 822.25 0 0 0 0 0 0 0 0 1 0 0
58 1164.17 880.00 0 0 0 0 0 0 0 0 0 1 0
59 968.67 729.50 0 0 0 0 0 0 0 0 0 0 1
60 888.61 778.00 0 0 0 0 0 0 0 0 0 0 0
> k <- length(x[1,])
> df <- as.data.frame(x)
> (mylm <- lm(df))
Call:
lm(formula = df)
Coefficients:
(Intercept) Gold M1 M2 M3 M4
966.1300 0.4304 79.2883 60.6833 44.0737 45.6598
M5 M6 M7 M8 M9 M10
63.6943 78.8141 48.7047 40.1346 56.9107 40.5558
M11
19.6090
> (mysum <- summary(mylm))
Call:
lm(formula = df)
Residuals:
Min 1Q Median 3Q Max
-412.34 -67.19 -10.31 68.88 219.96
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 966.1300 87.7929 11.005 1.33e-14 ***
Gold 0.4304 0.1006 4.278 9.18e-05 ***
M1 79.2883 85.7295 0.925 0.360
M2 60.6833 85.6080 0.709 0.482
M3 44.0737 85.5394 0.515 0.609
M4 45.6598 85.5857 0.533 0.596
M5 63.6943 85.5705 0.744 0.460
M6 78.8141 85.5817 0.921 0.362
M7 48.7047 85.5469 0.569 0.572
M8 40.1346 85.5492 0.469 0.641
M9 56.9107 85.6079 0.665 0.509
M10 40.5558 85.5198 0.474 0.638
M11 19.6090 85.5594 0.229 0.820
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 135.2 on 47 degrees of freedom
Multiple R-squared: 0.2908, Adjusted R-squared: 0.1097
F-statistic: 1.606 on 12 and 47 DF, p-value: 0.1226
> 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.614773e-02 5.229546e-02 0.9738523
[2,] 6.945027e-03 1.389005e-02 0.9930550
[3,] 1.897435e-03 3.794869e-03 0.9981026
[4,] 4.510338e-04 9.020677e-04 0.9995490
[5,] 1.499009e-04 2.998018e-04 0.9998501
[6,] 4.446916e-05 8.893833e-05 0.9999555
[7,] 8.664458e-06 1.732892e-05 0.9999913
[8,] 1.572618e-06 3.145237e-06 0.9999984
[9,] 3.687820e-07 7.375640e-07 0.9999996
[10,] 2.709001e-06 5.418001e-06 0.9999973
[11,] 2.206145e-06 4.412290e-06 0.9999978
[12,] 8.472691e-07 1.694538e-06 0.9999992
[13,] 2.011884e-07 4.023769e-07 0.9999998
[14,] 8.943963e-08 1.788793e-07 0.9999999
[15,] 8.034978e-08 1.606996e-07 0.9999999
[16,] 4.457841e-08 8.915682e-08 1.0000000
[17,] 1.993162e-08 3.986323e-08 1.0000000
[18,] 5.216782e-09 1.043356e-08 1.0000000
[19,] 2.107390e-09 4.214780e-09 1.0000000
[20,] 1.554537e-09 3.109074e-09 1.0000000
[21,] 7.459886e-10 1.491977e-09 1.0000000
[22,] 6.627961e-10 1.325592e-09 1.0000000
[23,] 3.615745e-10 7.231490e-10 1.0000000
[24,] 8.632745e-11 1.726549e-10 1.0000000
[25,] 3.295344e-11 6.590687e-11 1.0000000
[26,] 6.503253e-11 1.300651e-10 1.0000000
[27,] 2.239476e-10 4.478953e-10 1.0000000
[28,] 2.051352e-10 4.102704e-10 1.0000000
[29,] 4.291511e-11 8.583022e-11 1.0000000
> postscript(file="/var/www/html/rcomp/tmp/10mcg1259331118.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/2ap6a1259331118.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/3half1259331118.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/4sjve1259331118.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/53pju1259331118.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 = 60
Frequency = 1
1 2 3 4 5 6
-112.634730 -67.116328 -37.406780 -81.770184 -90.901816 -95.202180
7 8 9 10 11 12
-43.899921 -73.029723 -94.236048 -72.038428 -40.205537 12.762466
13 14 15 16 17 18
-17.583845 -26.680897 6.857750 -15.027148 -57.124431 -32.193175
19 20 21 22 23 24
-9.677448 42.151410 8.104018 21.534395 23.521584 68.278683
25 26 27 28 29 30
-25.217935 8.715874 27.626944 38.470150 -3.680846 -43.867697
31 32 33 34 35 36
-12.866064 -1.980164 13.485587 70.671265 127.738261 155.305364
37 38 39 40 41 42
97.290292 126.964470 108.084458 125.757154 162.656646 198.852452
43 44 45 46 47 48
208.048557 162.404638 161.718855 201.063486 219.960726 175.991875
49 50 51 52 53 54
58.146217 -41.883120 -105.162373 -67.429972 -10.949552 -27.589401
55 56 57 58 59 60
-141.605125 -129.546161 -89.072412 -221.230718 -331.015034 -412.338388
> postscript(file="/var/www/html/rcomp/tmp/642dp1259331118.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 = 60
Frequency = 1
lag(myerror, k = 1) myerror
0 -112.634730 NA
1 -67.116328 -112.634730
2 -37.406780 -67.116328
3 -81.770184 -37.406780
4 -90.901816 -81.770184
5 -95.202180 -90.901816
6 -43.899921 -95.202180
7 -73.029723 -43.899921
8 -94.236048 -73.029723
9 -72.038428 -94.236048
10 -40.205537 -72.038428
11 12.762466 -40.205537
12 -17.583845 12.762466
13 -26.680897 -17.583845
14 6.857750 -26.680897
15 -15.027148 6.857750
16 -57.124431 -15.027148
17 -32.193175 -57.124431
18 -9.677448 -32.193175
19 42.151410 -9.677448
20 8.104018 42.151410
21 21.534395 8.104018
22 23.521584 21.534395
23 68.278683 23.521584
24 -25.217935 68.278683
25 8.715874 -25.217935
26 27.626944 8.715874
27 38.470150 27.626944
28 -3.680846 38.470150
29 -43.867697 -3.680846
30 -12.866064 -43.867697
31 -1.980164 -12.866064
32 13.485587 -1.980164
33 70.671265 13.485587
34 127.738261 70.671265
35 155.305364 127.738261
36 97.290292 155.305364
37 126.964470 97.290292
38 108.084458 126.964470
39 125.757154 108.084458
40 162.656646 125.757154
41 198.852452 162.656646
42 208.048557 198.852452
43 162.404638 208.048557
44 161.718855 162.404638
45 201.063486 161.718855
46 219.960726 201.063486
47 175.991875 219.960726
48 58.146217 175.991875
49 -41.883120 58.146217
50 -105.162373 -41.883120
51 -67.429972 -105.162373
52 -10.949552 -67.429972
53 -27.589401 -10.949552
54 -141.605125 -27.589401
55 -129.546161 -141.605125
56 -89.072412 -129.546161
57 -221.230718 -89.072412
58 -331.015034 -221.230718
59 -412.338388 -331.015034
60 NA -412.338388
> dum1 <- dum[2:length(myerror),]
> dum1
lag(myerror, k = 1) myerror
[1,] -67.116328 -112.634730
[2,] -37.406780 -67.116328
[3,] -81.770184 -37.406780
[4,] -90.901816 -81.770184
[5,] -95.202180 -90.901816
[6,] -43.899921 -95.202180
[7,] -73.029723 -43.899921
[8,] -94.236048 -73.029723
[9,] -72.038428 -94.236048
[10,] -40.205537 -72.038428
[11,] 12.762466 -40.205537
[12,] -17.583845 12.762466
[13,] -26.680897 -17.583845
[14,] 6.857750 -26.680897
[15,] -15.027148 6.857750
[16,] -57.124431 -15.027148
[17,] -32.193175 -57.124431
[18,] -9.677448 -32.193175
[19,] 42.151410 -9.677448
[20,] 8.104018 42.151410
[21,] 21.534395 8.104018
[22,] 23.521584 21.534395
[23,] 68.278683 23.521584
[24,] -25.217935 68.278683
[25,] 8.715874 -25.217935
[26,] 27.626944 8.715874
[27,] 38.470150 27.626944
[28,] -3.680846 38.470150
[29,] -43.867697 -3.680846
[30,] -12.866064 -43.867697
[31,] -1.980164 -12.866064
[32,] 13.485587 -1.980164
[33,] 70.671265 13.485587
[34,] 127.738261 70.671265
[35,] 155.305364 127.738261
[36,] 97.290292 155.305364
[37,] 126.964470 97.290292
[38,] 108.084458 126.964470
[39,] 125.757154 108.084458
[40,] 162.656646 125.757154
[41,] 198.852452 162.656646
[42,] 208.048557 198.852452
[43,] 162.404638 208.048557
[44,] 161.718855 162.404638
[45,] 201.063486 161.718855
[46,] 219.960726 201.063486
[47,] 175.991875 219.960726
[48,] 58.146217 175.991875
[49,] -41.883120 58.146217
[50,] -105.162373 -41.883120
[51,] -67.429972 -105.162373
[52,] -10.949552 -67.429972
[53,] -27.589401 -10.949552
[54,] -141.605125 -27.589401
[55,] -129.546161 -141.605125
[56,] -89.072412 -129.546161
[57,] -221.230718 -89.072412
[58,] -331.015034 -221.230718
[59,] -412.338388 -331.015034
> z <- as.data.frame(dum1)
> z
lag(myerror, k = 1) myerror
1 -67.116328 -112.634730
2 -37.406780 -67.116328
3 -81.770184 -37.406780
4 -90.901816 -81.770184
5 -95.202180 -90.901816
6 -43.899921 -95.202180
7 -73.029723 -43.899921
8 -94.236048 -73.029723
9 -72.038428 -94.236048
10 -40.205537 -72.038428
11 12.762466 -40.205537
12 -17.583845 12.762466
13 -26.680897 -17.583845
14 6.857750 -26.680897
15 -15.027148 6.857750
16 -57.124431 -15.027148
17 -32.193175 -57.124431
18 -9.677448 -32.193175
19 42.151410 -9.677448
20 8.104018 42.151410
21 21.534395 8.104018
22 23.521584 21.534395
23 68.278683 23.521584
24 -25.217935 68.278683
25 8.715874 -25.217935
26 27.626944 8.715874
27 38.470150 27.626944
28 -3.680846 38.470150
29 -43.867697 -3.680846
30 -12.866064 -43.867697
31 -1.980164 -12.866064
32 13.485587 -1.980164
33 70.671265 13.485587
34 127.738261 70.671265
35 155.305364 127.738261
36 97.290292 155.305364
37 126.964470 97.290292
38 108.084458 126.964470
39 125.757154 108.084458
40 162.656646 125.757154
41 198.852452 162.656646
42 208.048557 198.852452
43 162.404638 208.048557
44 161.718855 162.404638
45 201.063486 161.718855
46 219.960726 201.063486
47 175.991875 219.960726
48 58.146217 175.991875
49 -41.883120 58.146217
50 -105.162373 -41.883120
51 -67.429972 -105.162373
52 -10.949552 -67.429972
53 -27.589401 -10.949552
54 -141.605125 -27.589401
55 -129.546161 -141.605125
56 -89.072412 -129.546161
57 -221.230718 -89.072412
58 -331.015034 -221.230718
59 -412.338388 -331.015034
> 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/7ehgi1259331118.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/8z0xd1259331118.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/9vljv1259331118.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/10lley1259331118.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/11xdcr1259331118.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/12l4gt1259331118.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/13tph81259331118.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/14ks0m1259331118.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/15odkl1259331118.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/1625w91259331118.tab")
+ }
>
> system("convert tmp/10mcg1259331118.ps tmp/10mcg1259331118.png")
> system("convert tmp/2ap6a1259331118.ps tmp/2ap6a1259331118.png")
> system("convert tmp/3half1259331118.ps tmp/3half1259331118.png")
> system("convert tmp/4sjve1259331118.ps tmp/4sjve1259331118.png")
> system("convert tmp/53pju1259331118.ps tmp/53pju1259331118.png")
> system("convert tmp/642dp1259331118.ps tmp/642dp1259331118.png")
> system("convert tmp/7ehgi1259331118.ps tmp/7ehgi1259331118.png")
> system("convert tmp/8z0xd1259331118.ps tmp/8z0xd1259331118.png")
> system("convert tmp/9vljv1259331118.ps tmp/9vljv1259331118.png")
> system("convert tmp/10lley1259331118.ps tmp/10lley1259331118.png")
>
>
> proc.time()
user system elapsed
2.498 1.624 5.312