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(10001.60,49.14,10411.75,44.61,10673.38,40.22,10539.51,44.23,10723.78,45.85,10682.06,53.38,10283.19,53.26,10377.18,51.8,10486.64,55.3,10545.38,57.81,10554.27,63.96,10532.54,63.77,10324.31,59.15,10695.25,56.12,10827.81,57.42,10872.48,63.52,10971.19,61.71,11145.65,63.01,11234.68,68.18,11333.88,72.03,10997.97,69.75,11036.89,74.41,11257.35,74.33,11533.59,64.24,11963.12,60.03,12185.15,59.44,12377.62,62.5,12512.89,55.04,12631.48,58.34,12268.53,61.92,12754.80,67.65,13407.75,67.68,13480.21,70.3,13673.28,75.26,13239.71,71.44,13557.69,76.36,13901.28,81.71,13200.58,92.6,13406.97,90.6,12538.12,92.23,12419.57,94.09,12193.88,102.79,12656.63,109.65,12812.48,124.05,12056.67,132.69,11322.38,135.81,11530.75,116.07,11114.08,101.42,9181.73,75.73,8614.55,55.48,8595.56,43.8,8396.20,45.29,7690.50,44.01,7235.47,47.48,7992.12,51.07,8398.37,57.84,8593.01,69.04,8679.75,65.61,9374.63,72.87,9634.97,68.41),dim=c(2,60),dimnames=list(c('Dow','Olie
'),1:60))
> y <- array(NA,dim=c(2,60),dimnames=list(c('Dow','Olie
'),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 = 'Linear Trend'
> par2 = 'Include Monthly Dummies'
> par1 = '1'
> #'GNU S' R Code compiled by R2WASP v. 1.0.44 ()
> #Author: Prof. Dr. P. Wessa
> #To cite this work: AUTHOR(S), (YEAR), YOUR SOFTWARE TITLE (vNUMBER) in Free Statistics Software (v$_version), Office for Research Development and Education, URL http://www.wessa.net/rwasp_YOURPAGE.wasp/
> #Source of accompanying publication: Office for Research, Development, and Education
> #Technical description: Write here your technical program description (don't use hard returns!)
> library(lattice)
> library(lmtest)
Loading required package: zoo
Attaching package: 'zoo'
The following object(s) are masked from package:base :
as.Date.numeric
> n25 <- 25 #minimum number of obs. for Goldfeld-Quandt test
> par1 <- as.numeric(par1)
> x <- t(y)
> k <- length(x[1,])
> n <- length(x[,1])
> x1 <- cbind(x[,par1], x[,1:k!=par1])
> mycolnames <- c(colnames(x)[par1], colnames(x)[1:k!=par1])
> colnames(x1) <- mycolnames #colnames(x)[par1]
> x <- x1
> if (par3 == 'First Differences'){
+ x2 <- array(0, dim=c(n-1,k), dimnames=list(1:(n-1), paste('(1-B)',colnames(x),sep='')))
+ for (i in 1:n-1) {
+ for (j in 1:k) {
+ x2[i,j] <- x[i+1,j] - x[i,j]
+ }
+ }
+ x <- x2
+ }
> if (par2 == 'Include Monthly Dummies'){
+ x2 <- array(0, dim=c(n,11), dimnames=list(1:n, paste('M', seq(1:11), sep ='')))
+ for (i in 1:11){
+ x2[seq(i,n,12),i] <- 1
+ }
+ x <- cbind(x, x2)
+ }
> if (par2 == 'Include Quarterly Dummies'){
+ x2 <- array(0, dim=c(n,3), dimnames=list(1:n, paste('Q', seq(1:3), sep ='')))
+ for (i in 1:3){
+ x2[seq(i,n,4),i] <- 1
+ }
+ x <- cbind(x, x2)
+ }
> k <- length(x[1,])
> if (par3 == 'Linear Trend'){
+ x <- cbind(x, c(1:n))
+ colnames(x)[k+1] <- 't'
+ }
> x
Dow Olie\r M1 M2 M3 M4 M5 M6 M7 M8 M9 M10 M11 t
1 10001.60 49.14 1 0 0 0 0 0 0 0 0 0 0 1
2 10411.75 44.61 0 1 0 0 0 0 0 0 0 0 0 2
3 10673.38 40.22 0 0 1 0 0 0 0 0 0 0 0 3
4 10539.51 44.23 0 0 0 1 0 0 0 0 0 0 0 4
5 10723.78 45.85 0 0 0 0 1 0 0 0 0 0 0 5
6 10682.06 53.38 0 0 0 0 0 1 0 0 0 0 0 6
7 10283.19 53.26 0 0 0 0 0 0 1 0 0 0 0 7
8 10377.18 51.80 0 0 0 0 0 0 0 1 0 0 0 8
9 10486.64 55.30 0 0 0 0 0 0 0 0 1 0 0 9
10 10545.38 57.81 0 0 0 0 0 0 0 0 0 1 0 10
11 10554.27 63.96 0 0 0 0 0 0 0 0 0 0 1 11
12 10532.54 63.77 0 0 0 0 0 0 0 0 0 0 0 12
13 10324.31 59.15 1 0 0 0 0 0 0 0 0 0 0 13
14 10695.25 56.12 0 1 0 0 0 0 0 0 0 0 0 14
15 10827.81 57.42 0 0 1 0 0 0 0 0 0 0 0 15
16 10872.48 63.52 0 0 0 1 0 0 0 0 0 0 0 16
17 10971.19 61.71 0 0 0 0 1 0 0 0 0 0 0 17
18 11145.65 63.01 0 0 0 0 0 1 0 0 0 0 0 18
19 11234.68 68.18 0 0 0 0 0 0 1 0 0 0 0 19
20 11333.88 72.03 0 0 0 0 0 0 0 1 0 0 0 20
21 10997.97 69.75 0 0 0 0 0 0 0 0 1 0 0 21
22 11036.89 74.41 0 0 0 0 0 0 0 0 0 1 0 22
23 11257.35 74.33 0 0 0 0 0 0 0 0 0 0 1 23
24 11533.59 64.24 0 0 0 0 0 0 0 0 0 0 0 24
25 11963.12 60.03 1 0 0 0 0 0 0 0 0 0 0 25
26 12185.15 59.44 0 1 0 0 0 0 0 0 0 0 0 26
27 12377.62 62.50 0 0 1 0 0 0 0 0 0 0 0 27
28 12512.89 55.04 0 0 0 1 0 0 0 0 0 0 0 28
29 12631.48 58.34 0 0 0 0 1 0 0 0 0 0 0 29
30 12268.53 61.92 0 0 0 0 0 1 0 0 0 0 0 30
31 12754.80 67.65 0 0 0 0 0 0 1 0 0 0 0 31
32 13407.75 67.68 0 0 0 0 0 0 0 1 0 0 0 32
33 13480.21 70.30 0 0 0 0 0 0 0 0 1 0 0 33
34 13673.28 75.26 0 0 0 0 0 0 0 0 0 1 0 34
35 13239.71 71.44 0 0 0 0 0 0 0 0 0 0 1 35
36 13557.69 76.36 0 0 0 0 0 0 0 0 0 0 0 36
37 13901.28 81.71 1 0 0 0 0 0 0 0 0 0 0 37
38 13200.58 92.60 0 1 0 0 0 0 0 0 0 0 0 38
39 13406.97 90.60 0 0 1 0 0 0 0 0 0 0 0 39
40 12538.12 92.23 0 0 0 1 0 0 0 0 0 0 0 40
41 12419.57 94.09 0 0 0 0 1 0 0 0 0 0 0 41
42 12193.88 102.79 0 0 0 0 0 1 0 0 0 0 0 42
43 12656.63 109.65 0 0 0 0 0 0 1 0 0 0 0 43
44 12812.48 124.05 0 0 0 0 0 0 0 1 0 0 0 44
45 12056.67 132.69 0 0 0 0 0 0 0 0 1 0 0 45
46 11322.38 135.81 0 0 0 0 0 0 0 0 0 1 0 46
47 11530.75 116.07 0 0 0 0 0 0 0 0 0 0 1 47
48 11114.08 101.42 0 0 0 0 0 0 0 0 0 0 0 48
49 9181.73 75.73 1 0 0 0 0 0 0 0 0 0 0 49
50 8614.55 55.48 0 1 0 0 0 0 0 0 0 0 0 50
51 8595.56 43.80 0 0 1 0 0 0 0 0 0 0 0 51
52 8396.20 45.29 0 0 0 1 0 0 0 0 0 0 0 52
53 7690.50 44.01 0 0 0 0 1 0 0 0 0 0 0 53
54 7235.47 47.48 0 0 0 0 0 1 0 0 0 0 0 54
55 7992.12 51.07 0 0 0 0 0 0 1 0 0 0 0 55
56 8398.37 57.84 0 0 0 0 0 0 0 1 0 0 0 56
57 8593.01 69.04 0 0 0 0 0 0 0 0 1 0 0 57
58 8679.75 65.61 0 0 0 0 0 0 0 0 0 1 0 58
59 9374.63 72.87 0 0 0 0 0 0 0 0 0 0 1 59
60 9634.97 68.41 0 0 0 0 0 0 0 0 0 0 0 60
> k <- length(x[1,])
> df <- as.data.frame(x)
> (mylm <- lm(df))
Call:
lm(formula = df)
Coefficients:
(Intercept) `Olie\r` M1 M2 M3 M4
8856.12 53.74 -169.42 10.34 357.03 135.12
M5 M6 M7 M8 M9 M10
55.46 -346.36 -250.83 -178.17 -531.17 -685.03
M11 t
-390.75 -44.53
> (mysum <- summary(mylm))
Call:
lm(formula = df)
Residuals:
Min 1Q Median 3Q Max
-2098.1 -861.7 -419.0 531.1 2972.1
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 8856.116 927.094 9.553 1.71e-12 ***
`Olie\r` 53.736 9.631 5.579 1.23e-06 ***
M1 -169.424 893.017 -0.190 0.850363
M2 10.341 894.472 0.012 0.990826
M3 357.029 896.690 0.398 0.692351
M4 135.122 894.761 0.151 0.880625
M5 55.461 893.472 0.062 0.950773
M6 -346.360 888.590 -0.390 0.698494
M7 -250.826 886.072 -0.283 0.778389
M8 -178.173 885.260 -0.201 0.841378
M9 -531.167 886.502 -0.599 0.551998
M10 -685.032 887.562 -0.772 0.444173
M11 -390.749 885.658 -0.441 0.661138
t -44.532 11.454 -3.888 0.000323 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 1398 on 46 degrees of freedom
Multiple R-squared: 0.437, Adjusted R-squared: 0.2779
F-statistic: 2.746 on 13 and 46 DF, p-value: 0.005827
> 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,] 1.032106e-04 2.064213e-04 0.9998968
[2,] 1.329817e-05 2.659635e-05 0.9999867
[3,] 2.654261e-04 5.308522e-04 0.9997346
[4,] 1.702104e-04 3.404208e-04 0.9998298
[5,] 3.910236e-05 7.820471e-05 0.9999609
[6,] 1.059097e-05 2.118195e-05 0.9999894
[7,] 8.144205e-06 1.628841e-05 0.9999919
[8,] 1.874737e-05 3.749473e-05 0.9999813
[9,] 5.611520e-05 1.122304e-04 0.9999439
[10,] 3.747949e-05 7.495898e-05 0.9999625
[11,] 6.045054e-05 1.209011e-04 0.9999395
[12,] 2.569432e-05 5.138865e-05 0.9999743
[13,] 9.238759e-06 1.847752e-05 0.9999908
[14,] 5.841993e-06 1.168399e-05 0.9999942
[15,] 6.956817e-06 1.391363e-05 0.9999930
[16,] 4.199038e-05 8.398076e-05 0.9999580
[17,] 2.281943e-04 4.563885e-04 0.9997718
[18,] 9.344578e-04 1.868916e-03 0.9990655
[19,] 4.633461e-04 9.266922e-04 0.9995367
[20,] 1.569279e-03 3.138558e-03 0.9984307
[21,] 6.348968e-03 1.269794e-02 0.9936510
[22,] 6.677943e-03 1.335589e-02 0.9933221
[23,] 8.145583e-03 1.629117e-02 0.9918544
[24,] 7.150239e-03 1.430048e-02 0.9928498
[25,] 1.645192e-02 3.290384e-02 0.9835481
[26,] 6.299726e-02 1.259945e-01 0.9370027
[27,] 1.988311e-01 3.976622e-01 0.8011689
> postscript(file="/var/www/html/rcomp/tmp/1e14c1262113898.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/2ryeu1262113898.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/3zy541262113898.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/4bow51262113898.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/5p1iz1262113898.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
-1281.15674 -762.81392 -567.43866 -650.35126 -428.94059 -428.94102
7 8 9 10 11 12
-872.36476 -728.04050 -409.13056 -286.87204 -858.20996 -1215.94705
13 14 15 16 17 18
-961.95972 -563.43121 -802.88494 -819.56619 -499.40035 51.55576
19 20 21 22 23 24
-188.23249 -324.03747 -139.90227 -152.99659 -177.98797 294.23334
25 26 27 28 29 30
1163.94883 1282.45101 1008.33156 1810.91321 1876.36706 1767.39462
31 32 33 34 35 36
1894.75410 2517.97141 2847.16922 2972.10404 2494.05606 2201.43697
37 38 39 40 41 42
2471.49437 1050.37495 1062.08069 372.08025 277.77424 30.93245
43 44 45 46 47 48
74.05001 -572.02189 -1394.58603 -2098.13660 -1078.76424 -1054.41585
49 50 51 52 53 54
-1392.32674 -1006.58082 -700.08865 -713.07601 -1225.80036 -1420.94181
55 56 57 58 59 60
-908.20686 -893.87155 -903.55036 -434.09880 -379.09390 -225.30741
> postscript(file="/var/www/html/rcomp/tmp/6jl8f1262113898.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 -1281.15674 NA
1 -762.81392 -1281.15674
2 -567.43866 -762.81392
3 -650.35126 -567.43866
4 -428.94059 -650.35126
5 -428.94102 -428.94059
6 -872.36476 -428.94102
7 -728.04050 -872.36476
8 -409.13056 -728.04050
9 -286.87204 -409.13056
10 -858.20996 -286.87204
11 -1215.94705 -858.20996
12 -961.95972 -1215.94705
13 -563.43121 -961.95972
14 -802.88494 -563.43121
15 -819.56619 -802.88494
16 -499.40035 -819.56619
17 51.55576 -499.40035
18 -188.23249 51.55576
19 -324.03747 -188.23249
20 -139.90227 -324.03747
21 -152.99659 -139.90227
22 -177.98797 -152.99659
23 294.23334 -177.98797
24 1163.94883 294.23334
25 1282.45101 1163.94883
26 1008.33156 1282.45101
27 1810.91321 1008.33156
28 1876.36706 1810.91321
29 1767.39462 1876.36706
30 1894.75410 1767.39462
31 2517.97141 1894.75410
32 2847.16922 2517.97141
33 2972.10404 2847.16922
34 2494.05606 2972.10404
35 2201.43697 2494.05606
36 2471.49437 2201.43697
37 1050.37495 2471.49437
38 1062.08069 1050.37495
39 372.08025 1062.08069
40 277.77424 372.08025
41 30.93245 277.77424
42 74.05001 30.93245
43 -572.02189 74.05001
44 -1394.58603 -572.02189
45 -2098.13660 -1394.58603
46 -1078.76424 -2098.13660
47 -1054.41585 -1078.76424
48 -1392.32674 -1054.41585
49 -1006.58082 -1392.32674
50 -700.08865 -1006.58082
51 -713.07601 -700.08865
52 -1225.80036 -713.07601
53 -1420.94181 -1225.80036
54 -908.20686 -1420.94181
55 -893.87155 -908.20686
56 -903.55036 -893.87155
57 -434.09880 -903.55036
58 -379.09390 -434.09880
59 -225.30741 -379.09390
60 NA -225.30741
> dum1 <- dum[2:length(myerror),]
> dum1
lag(myerror, k = 1) myerror
[1,] -762.81392 -1281.15674
[2,] -567.43866 -762.81392
[3,] -650.35126 -567.43866
[4,] -428.94059 -650.35126
[5,] -428.94102 -428.94059
[6,] -872.36476 -428.94102
[7,] -728.04050 -872.36476
[8,] -409.13056 -728.04050
[9,] -286.87204 -409.13056
[10,] -858.20996 -286.87204
[11,] -1215.94705 -858.20996
[12,] -961.95972 -1215.94705
[13,] -563.43121 -961.95972
[14,] -802.88494 -563.43121
[15,] -819.56619 -802.88494
[16,] -499.40035 -819.56619
[17,] 51.55576 -499.40035
[18,] -188.23249 51.55576
[19,] -324.03747 -188.23249
[20,] -139.90227 -324.03747
[21,] -152.99659 -139.90227
[22,] -177.98797 -152.99659
[23,] 294.23334 -177.98797
[24,] 1163.94883 294.23334
[25,] 1282.45101 1163.94883
[26,] 1008.33156 1282.45101
[27,] 1810.91321 1008.33156
[28,] 1876.36706 1810.91321
[29,] 1767.39462 1876.36706
[30,] 1894.75410 1767.39462
[31,] 2517.97141 1894.75410
[32,] 2847.16922 2517.97141
[33,] 2972.10404 2847.16922
[34,] 2494.05606 2972.10404
[35,] 2201.43697 2494.05606
[36,] 2471.49437 2201.43697
[37,] 1050.37495 2471.49437
[38,] 1062.08069 1050.37495
[39,] 372.08025 1062.08069
[40,] 277.77424 372.08025
[41,] 30.93245 277.77424
[42,] 74.05001 30.93245
[43,] -572.02189 74.05001
[44,] -1394.58603 -572.02189
[45,] -2098.13660 -1394.58603
[46,] -1078.76424 -2098.13660
[47,] -1054.41585 -1078.76424
[48,] -1392.32674 -1054.41585
[49,] -1006.58082 -1392.32674
[50,] -700.08865 -1006.58082
[51,] -713.07601 -700.08865
[52,] -1225.80036 -713.07601
[53,] -1420.94181 -1225.80036
[54,] -908.20686 -1420.94181
[55,] -893.87155 -908.20686
[56,] -903.55036 -893.87155
[57,] -434.09880 -903.55036
[58,] -379.09390 -434.09880
[59,] -225.30741 -379.09390
> z <- as.data.frame(dum1)
> z
lag(myerror, k = 1) myerror
1 -762.81392 -1281.15674
2 -567.43866 -762.81392
3 -650.35126 -567.43866
4 -428.94059 -650.35126
5 -428.94102 -428.94059
6 -872.36476 -428.94102
7 -728.04050 -872.36476
8 -409.13056 -728.04050
9 -286.87204 -409.13056
10 -858.20996 -286.87204
11 -1215.94705 -858.20996
12 -961.95972 -1215.94705
13 -563.43121 -961.95972
14 -802.88494 -563.43121
15 -819.56619 -802.88494
16 -499.40035 -819.56619
17 51.55576 -499.40035
18 -188.23249 51.55576
19 -324.03747 -188.23249
20 -139.90227 -324.03747
21 -152.99659 -139.90227
22 -177.98797 -152.99659
23 294.23334 -177.98797
24 1163.94883 294.23334
25 1282.45101 1163.94883
26 1008.33156 1282.45101
27 1810.91321 1008.33156
28 1876.36706 1810.91321
29 1767.39462 1876.36706
30 1894.75410 1767.39462
31 2517.97141 1894.75410
32 2847.16922 2517.97141
33 2972.10404 2847.16922
34 2494.05606 2972.10404
35 2201.43697 2494.05606
36 2471.49437 2201.43697
37 1050.37495 2471.49437
38 1062.08069 1050.37495
39 372.08025 1062.08069
40 277.77424 372.08025
41 30.93245 277.77424
42 74.05001 30.93245
43 -572.02189 74.05001
44 -1394.58603 -572.02189
45 -2098.13660 -1394.58603
46 -1078.76424 -2098.13660
47 -1054.41585 -1078.76424
48 -1392.32674 -1054.41585
49 -1006.58082 -1392.32674
50 -700.08865 -1006.58082
51 -713.07601 -700.08865
52 -1225.80036 -713.07601
53 -1420.94181 -1225.80036
54 -908.20686 -1420.94181
55 -893.87155 -908.20686
56 -903.55036 -893.87155
57 -434.09880 -903.55036
58 -379.09390 -434.09880
59 -225.30741 -379.09390
> 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/7e67p1262113898.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/8mozy1262113898.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/9l0hx1262113898.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/10nwea1262113898.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/11gkgv1262113898.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/128a8b1262113898.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/13313h1262113898.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/14wnfs1262113898.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/15vgnz1262113898.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/16xlf91262113898.tab")
+ }
>
> try(system("convert tmp/1e14c1262113898.ps tmp/1e14c1262113898.png",intern=TRUE))
character(0)
> try(system("convert tmp/2ryeu1262113898.ps tmp/2ryeu1262113898.png",intern=TRUE))
character(0)
> try(system("convert tmp/3zy541262113898.ps tmp/3zy541262113898.png",intern=TRUE))
character(0)
> try(system("convert tmp/4bow51262113898.ps tmp/4bow51262113898.png",intern=TRUE))
character(0)
> try(system("convert tmp/5p1iz1262113898.ps tmp/5p1iz1262113898.png",intern=TRUE))
character(0)
> try(system("convert tmp/6jl8f1262113898.ps tmp/6jl8f1262113898.png",intern=TRUE))
character(0)
> try(system("convert tmp/7e67p1262113898.ps tmp/7e67p1262113898.png",intern=TRUE))
character(0)
> try(system("convert tmp/8mozy1262113898.ps tmp/8mozy1262113898.png",intern=TRUE))
character(0)
> try(system("convert tmp/9l0hx1262113898.ps tmp/9l0hx1262113898.png",intern=TRUE))
character(0)
> try(system("convert tmp/10nwea1262113898.ps tmp/10nwea1262113898.png",intern=TRUE))
character(0)
>
>
> proc.time()
user system elapsed
2.364 1.533 2.924