R version 2.8.0 (2008-10-20)
Copyright (C) 2008 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.
Natural language support but running in an English locale
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(10554.27,10532.54,10324.31,10695.25,10827.81,10872.48,10971.19,11145.65,11234.68,11333.88,10997.97,11036.89,11257.35,11533.59,11963.12,12185.15,12377.62,12512.89,12631.48,12268.53,12754.8,13407.75,13480.21,13673.28,13239.71,13557.69,13901.28,13200.58,13406.97,12538.12,12419.57,12193.88,12656.63,12812.48,12056.67,11322.38,11530.75,11114.08,9181.73,8614.55,8595.56,8396.2,7690.5,7235.47,7992.12,8398.37,8593,8679.75,9374.63,9634.97,9857.34,10238.83,10433.44,10471.24,10214.51,10677.52,11052.15,10500.19,10159.27,10222.24,10350.4),dim=c(1,61),dimnames=list(c('dowjones'),1:61))
> y <- array(NA,dim=c(1,61),dimnames=list(c('dowjones'),1:61))
> 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
dowjones M1 M2 M3 M4 M5 M6 M7 M8 M9 M10 M11 t
1 10554.27 1 0 0 0 0 0 0 0 0 0 0 1
2 10532.54 0 1 0 0 0 0 0 0 0 0 0 2
3 10324.31 0 0 1 0 0 0 0 0 0 0 0 3
4 10695.25 0 0 0 1 0 0 0 0 0 0 0 4
5 10827.81 0 0 0 0 1 0 0 0 0 0 0 5
6 10872.48 0 0 0 0 0 1 0 0 0 0 0 6
7 10971.19 0 0 0 0 0 0 1 0 0 0 0 7
8 11145.65 0 0 0 0 0 0 0 1 0 0 0 8
9 11234.68 0 0 0 0 0 0 0 0 1 0 0 9
10 11333.88 0 0 0 0 0 0 0 0 0 1 0 10
11 10997.97 0 0 0 0 0 0 0 0 0 0 1 11
12 11036.89 0 0 0 0 0 0 0 0 0 0 0 12
13 11257.35 1 0 0 0 0 0 0 0 0 0 0 13
14 11533.59 0 1 0 0 0 0 0 0 0 0 0 14
15 11963.12 0 0 1 0 0 0 0 0 0 0 0 15
16 12185.15 0 0 0 1 0 0 0 0 0 0 0 16
17 12377.62 0 0 0 0 1 0 0 0 0 0 0 17
18 12512.89 0 0 0 0 0 1 0 0 0 0 0 18
19 12631.48 0 0 0 0 0 0 1 0 0 0 0 19
20 12268.53 0 0 0 0 0 0 0 1 0 0 0 20
21 12754.80 0 0 0 0 0 0 0 0 1 0 0 21
22 13407.75 0 0 0 0 0 0 0 0 0 1 0 22
23 13480.21 0 0 0 0 0 0 0 0 0 0 1 23
24 13673.28 0 0 0 0 0 0 0 0 0 0 0 24
25 13239.71 1 0 0 0 0 0 0 0 0 0 0 25
26 13557.69 0 1 0 0 0 0 0 0 0 0 0 26
27 13901.28 0 0 1 0 0 0 0 0 0 0 0 27
28 13200.58 0 0 0 1 0 0 0 0 0 0 0 28
29 13406.97 0 0 0 0 1 0 0 0 0 0 0 29
30 12538.12 0 0 0 0 0 1 0 0 0 0 0 30
31 12419.57 0 0 0 0 0 0 1 0 0 0 0 31
32 12193.88 0 0 0 0 0 0 0 1 0 0 0 32
33 12656.63 0 0 0 0 0 0 0 0 1 0 0 33
34 12812.48 0 0 0 0 0 0 0 0 0 1 0 34
35 12056.67 0 0 0 0 0 0 0 0 0 0 1 35
36 11322.38 0 0 0 0 0 0 0 0 0 0 0 36
37 11530.75 1 0 0 0 0 0 0 0 0 0 0 37
38 11114.08 0 1 0 0 0 0 0 0 0 0 0 38
39 9181.73 0 0 1 0 0 0 0 0 0 0 0 39
40 8614.55 0 0 0 1 0 0 0 0 0 0 0 40
41 8595.56 0 0 0 0 1 0 0 0 0 0 0 41
42 8396.20 0 0 0 0 0 1 0 0 0 0 0 42
43 7690.50 0 0 0 0 0 0 1 0 0 0 0 43
44 7235.47 0 0 0 0 0 0 0 1 0 0 0 44
45 7992.12 0 0 0 0 0 0 0 0 1 0 0 45
46 8398.37 0 0 0 0 0 0 0 0 0 1 0 46
47 8593.00 0 0 0 0 0 0 0 0 0 0 1 47
48 8679.75 0 0 0 0 0 0 0 0 0 0 0 48
49 9374.63 1 0 0 0 0 0 0 0 0 0 0 49
50 9634.97 0 1 0 0 0 0 0 0 0 0 0 50
51 9857.34 0 0 1 0 0 0 0 0 0 0 0 51
52 10238.83 0 0 0 1 0 0 0 0 0 0 0 52
53 10433.44 0 0 0 0 1 0 0 0 0 0 0 53
54 10471.24 0 0 0 0 0 1 0 0 0 0 0 54
55 10214.51 0 0 0 0 0 0 1 0 0 0 0 55
56 10677.52 0 0 0 0 0 0 0 1 0 0 0 56
57 11052.15 0 0 0 0 0 0 0 0 1 0 0 57
58 10500.19 0 0 0 0 0 0 0 0 0 1 0 58
59 10159.27 0 0 0 0 0 0 0 0 0 0 1 59
60 10222.24 0 0 0 0 0 0 0 0 0 0 0 60
61 10350.40 1 0 0 0 0 0 0 0 0 0 0 61
> k <- length(x[1,])
> df <- as.data.frame(x)
> (mylm <- lm(df))
Call:
lm(formula = df)
Coefficients:
(Intercept) M1 M2 M3 M4 M5
12434.14 -136.73 -114.34 -303.16 -321.64 -140.03
M6 M7 M8 M9 M10 M11
-269.93 -402.46 -443.50 30.57 223.22 30.32
t
-40.20
> (mysum <- summary(mylm))
Call:
lm(formula = df)
Residuals:
Min 1Q Median 3Q Max
-2986.3 -1050.5 216.8 1081.9 2855.7
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 12434.14 857.66 14.498 < 2e-16 ***
M1 -136.73 1000.23 -0.137 0.89184
M2 -114.34 1049.85 -0.109 0.91373
M3 -303.16 1048.51 -0.289 0.77372
M4 -321.64 1047.30 -0.307 0.76008
M5 -140.03 1046.24 -0.134 0.89409
M6 -269.93 1045.32 -0.258 0.79734
M7 -402.46 1044.54 -0.385 0.70172
M8 -443.50 1043.91 -0.425 0.67285
M9 30.57 1043.41 0.029 0.97675
M10 223.22 1043.05 0.214 0.83145
M11 30.32 1042.84 0.029 0.97693
t -40.20 12.17 -3.304 0.00181 **
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Residual standard error: 1649 on 48 degrees of freedom
Multiple R-squared: 0.1935, Adjusted R-squared: -0.008181
F-statistic: 0.9594 on 12 and 48 DF, p-value: 0.499
> 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.001609e-02 2.003218e-02 0.989983909
[2,] 2.112749e-03 4.225498e-03 0.997887251
[3,] 4.615126e-04 9.230253e-04 0.999538487
[4,] 9.286127e-05 1.857225e-04 0.999907139
[5,] 1.649521e-05 3.299042e-05 0.999983505
[6,] 2.481545e-06 4.963090e-06 0.999997518
[7,] 1.736291e-06 3.472582e-06 0.999998264
[8,] 3.412654e-06 6.825308e-06 0.999996587
[9,] 4.787701e-06 9.575402e-06 0.999995212
[10,] 9.931454e-07 1.986291e-06 0.999999007
[11,] 2.171464e-07 4.342928e-07 0.999999783
[12,] 7.981900e-08 1.596380e-07 0.999999920
[13,] 6.433409e-08 1.286682e-07 0.999999936
[14,] 4.351694e-08 8.703387e-08 0.999999956
[15,] 4.128092e-07 8.256183e-07 0.999999587
[16,] 2.857600e-06 5.715199e-06 0.999997142
[17,] 1.278847e-05 2.557694e-05 0.999987212
[18,] 3.270886e-05 6.541772e-05 0.999967291
[19,] 2.015060e-04 4.030120e-04 0.999798494
[20,] 2.616149e-03 5.232298e-03 0.997383851
[21,] 4.404936e-02 8.809871e-02 0.955950645
[22,] 4.269181e-01 8.538362e-01 0.573081886
[23,] 9.400994e-01 1.198011e-01 0.059900563
[24,] 9.904007e-01 1.919866e-02 0.009599329
[25,] 9.933778e-01 1.324432e-02 0.006622161
[26,] 9.918264e-01 1.634727e-02 0.008173634
[27,] 9.854311e-01 2.913785e-02 0.014568926
[28,] 9.746186e-01 5.076281e-02 0.025381403
[29,] 9.826679e-01 3.466415e-02 0.017332077
[30,] 9.924273e-01 1.514543e-02 0.007572717
> postscript(file="/var/www/html/freestat/rcomp/tmp/16u1q1293385355.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/freestat/rcomp/tmp/26u1q1293385355.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/freestat/rcomp/tmp/3zl1c1293385355.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/freestat/rcomp/tmp/4zl1c1293385355.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/freestat/rcomp/tmp/5zl1c1293385355.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 = 61
Frequency = 1
1 2 3 4 5 6
-1702.93902 -1706.85322 -1686.06522 -1256.44122 -1265.28922 -1050.52522
7 8 9 10 11 12
-779.07922 -523.37922 -868.21522 -921.47322 -1024.27322 -914.83722
13 14 15 16 17 18
-517.44941 -223.39361 435.15439 715.86839 766.93039 1072.29439
19 20 21 22 23 24
1363.62039 1081.91039 1134.31439 1634.80639 1940.37639 2203.96239
25 26 27 28 29 30
1947.32020 2283.11600 2855.72400 2213.70800 2278.69000 1579.93400
31 32 33 34 35 36
1634.12000 1489.67000 1518.55400 1521.94600 999.24600 335.47200
37 38 39 40 41 42
720.76980 321.91561 -1381.41639 -1889.91239 -2050.31039 -2079.57639
43 44 45 46 47 48
-2612.54039 -2986.33039 -2663.54639 -2409.75439 -1982.01439 -1824.74839
49 50 51 52 53 54
-952.94059 -674.78478 -223.39678 216.77722 269.97922 477.87322
55 56 57 58 59 60
393.87922 938.12922 878.89322 174.47522 66.66522 200.15122
61
505.23902
> postscript(file="/var/www/html/freestat/rcomp/tmp/6acie1293385355.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 = 61
Frequency = 1
lag(myerror, k = 1) myerror
0 -1702.93902 NA
1 -1706.85322 -1702.93902
2 -1686.06522 -1706.85322
3 -1256.44122 -1686.06522
4 -1265.28922 -1256.44122
5 -1050.52522 -1265.28922
6 -779.07922 -1050.52522
7 -523.37922 -779.07922
8 -868.21522 -523.37922
9 -921.47322 -868.21522
10 -1024.27322 -921.47322
11 -914.83722 -1024.27322
12 -517.44941 -914.83722
13 -223.39361 -517.44941
14 435.15439 -223.39361
15 715.86839 435.15439
16 766.93039 715.86839
17 1072.29439 766.93039
18 1363.62039 1072.29439
19 1081.91039 1363.62039
20 1134.31439 1081.91039
21 1634.80639 1134.31439
22 1940.37639 1634.80639
23 2203.96239 1940.37639
24 1947.32020 2203.96239
25 2283.11600 1947.32020
26 2855.72400 2283.11600
27 2213.70800 2855.72400
28 2278.69000 2213.70800
29 1579.93400 2278.69000
30 1634.12000 1579.93400
31 1489.67000 1634.12000
32 1518.55400 1489.67000
33 1521.94600 1518.55400
34 999.24600 1521.94600
35 335.47200 999.24600
36 720.76980 335.47200
37 321.91561 720.76980
38 -1381.41639 321.91561
39 -1889.91239 -1381.41639
40 -2050.31039 -1889.91239
41 -2079.57639 -2050.31039
42 -2612.54039 -2079.57639
43 -2986.33039 -2612.54039
44 -2663.54639 -2986.33039
45 -2409.75439 -2663.54639
46 -1982.01439 -2409.75439
47 -1824.74839 -1982.01439
48 -952.94059 -1824.74839
49 -674.78478 -952.94059
50 -223.39678 -674.78478
51 216.77722 -223.39678
52 269.97922 216.77722
53 477.87322 269.97922
54 393.87922 477.87322
55 938.12922 393.87922
56 878.89322 938.12922
57 174.47522 878.89322
58 66.66522 174.47522
59 200.15122 66.66522
60 505.23902 200.15122
61 NA 505.23902
> dum1 <- dum[2:length(myerror),]
> dum1
lag(myerror, k = 1) myerror
[1,] -1706.85322 -1702.93902
[2,] -1686.06522 -1706.85322
[3,] -1256.44122 -1686.06522
[4,] -1265.28922 -1256.44122
[5,] -1050.52522 -1265.28922
[6,] -779.07922 -1050.52522
[7,] -523.37922 -779.07922
[8,] -868.21522 -523.37922
[9,] -921.47322 -868.21522
[10,] -1024.27322 -921.47322
[11,] -914.83722 -1024.27322
[12,] -517.44941 -914.83722
[13,] -223.39361 -517.44941
[14,] 435.15439 -223.39361
[15,] 715.86839 435.15439
[16,] 766.93039 715.86839
[17,] 1072.29439 766.93039
[18,] 1363.62039 1072.29439
[19,] 1081.91039 1363.62039
[20,] 1134.31439 1081.91039
[21,] 1634.80639 1134.31439
[22,] 1940.37639 1634.80639
[23,] 2203.96239 1940.37639
[24,] 1947.32020 2203.96239
[25,] 2283.11600 1947.32020
[26,] 2855.72400 2283.11600
[27,] 2213.70800 2855.72400
[28,] 2278.69000 2213.70800
[29,] 1579.93400 2278.69000
[30,] 1634.12000 1579.93400
[31,] 1489.67000 1634.12000
[32,] 1518.55400 1489.67000
[33,] 1521.94600 1518.55400
[34,] 999.24600 1521.94600
[35,] 335.47200 999.24600
[36,] 720.76980 335.47200
[37,] 321.91561 720.76980
[38,] -1381.41639 321.91561
[39,] -1889.91239 -1381.41639
[40,] -2050.31039 -1889.91239
[41,] -2079.57639 -2050.31039
[42,] -2612.54039 -2079.57639
[43,] -2986.33039 -2612.54039
[44,] -2663.54639 -2986.33039
[45,] -2409.75439 -2663.54639
[46,] -1982.01439 -2409.75439
[47,] -1824.74839 -1982.01439
[48,] -952.94059 -1824.74839
[49,] -674.78478 -952.94059
[50,] -223.39678 -674.78478
[51,] 216.77722 -223.39678
[52,] 269.97922 216.77722
[53,] 477.87322 269.97922
[54,] 393.87922 477.87322
[55,] 938.12922 393.87922
[56,] 878.89322 938.12922
[57,] 174.47522 878.89322
[58,] 66.66522 174.47522
[59,] 200.15122 66.66522
[60,] 505.23902 200.15122
> z <- as.data.frame(dum1)
> z
lag(myerror, k = 1) myerror
1 -1706.85322 -1702.93902
2 -1686.06522 -1706.85322
3 -1256.44122 -1686.06522
4 -1265.28922 -1256.44122
5 -1050.52522 -1265.28922
6 -779.07922 -1050.52522
7 -523.37922 -779.07922
8 -868.21522 -523.37922
9 -921.47322 -868.21522
10 -1024.27322 -921.47322
11 -914.83722 -1024.27322
12 -517.44941 -914.83722
13 -223.39361 -517.44941
14 435.15439 -223.39361
15 715.86839 435.15439
16 766.93039 715.86839
17 1072.29439 766.93039
18 1363.62039 1072.29439
19 1081.91039 1363.62039
20 1134.31439 1081.91039
21 1634.80639 1134.31439
22 1940.37639 1634.80639
23 2203.96239 1940.37639
24 1947.32020 2203.96239
25 2283.11600 1947.32020
26 2855.72400 2283.11600
27 2213.70800 2855.72400
28 2278.69000 2213.70800
29 1579.93400 2278.69000
30 1634.12000 1579.93400
31 1489.67000 1634.12000
32 1518.55400 1489.67000
33 1521.94600 1518.55400
34 999.24600 1521.94600
35 335.47200 999.24600
36 720.76980 335.47200
37 321.91561 720.76980
38 -1381.41639 321.91561
39 -1889.91239 -1381.41639
40 -2050.31039 -1889.91239
41 -2079.57639 -2050.31039
42 -2612.54039 -2079.57639
43 -2986.33039 -2612.54039
44 -2663.54639 -2986.33039
45 -2409.75439 -2663.54639
46 -1982.01439 -2409.75439
47 -1824.74839 -1982.01439
48 -952.94059 -1824.74839
49 -674.78478 -952.94059
50 -223.39678 -674.78478
51 216.77722 -223.39678
52 269.97922 216.77722
53 477.87322 269.97922
54 393.87922 477.87322
55 938.12922 393.87922
56 878.89322 938.12922
57 174.47522 878.89322
58 66.66522 174.47522
59 200.15122 66.66522
60 505.23902 200.15122
> 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/freestat/rcomp/tmp/7lmzz1293385355.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/freestat/rcomp/tmp/8lmzz1293385355.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/freestat/rcomp/tmp/9lmzz1293385355.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/freestat/rcomp/tmp/10vdgk1293385355.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/freestat/rcomp/createtable file can be downloaded at http://www.wessa.net/cretab
> load(file="/var/www/html/freestat/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/freestat/rcomp/tmp/11zdxq1293385355.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/freestat/rcomp/tmp/122www1293385355.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/freestat/rcomp/tmp/13y6tn1293385355.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/freestat/rcomp/tmp/14k6at1293385355.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/freestat/rcomp/tmp/15npqz1293385355.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/freestat/rcomp/tmp/1687p41293385355.tab")
+ }
> try(system("convert tmp/16u1q1293385355.ps tmp/16u1q1293385355.png",intern=TRUE))
character(0)
> try(system("convert tmp/26u1q1293385355.ps tmp/26u1q1293385355.png",intern=TRUE))
character(0)
> try(system("convert tmp/3zl1c1293385355.ps tmp/3zl1c1293385355.png",intern=TRUE))
character(0)
> try(system("convert tmp/4zl1c1293385355.ps tmp/4zl1c1293385355.png",intern=TRUE))
character(0)
> try(system("convert tmp/5zl1c1293385355.ps tmp/5zl1c1293385355.png",intern=TRUE))
character(0)
> try(system("convert tmp/6acie1293385355.ps tmp/6acie1293385355.png",intern=TRUE))
character(0)
> try(system("convert tmp/7lmzz1293385355.ps tmp/7lmzz1293385355.png",intern=TRUE))
character(0)
> try(system("convert tmp/8lmzz1293385355.ps tmp/8lmzz1293385355.png",intern=TRUE))
character(0)
> try(system("convert tmp/9lmzz1293385355.ps tmp/9lmzz1293385355.png",intern=TRUE))
character(0)
> try(system("convert tmp/10vdgk1293385355.ps tmp/10vdgk1293385355.png",intern=TRUE))
character(0)
>
>
> proc.time()
user system elapsed
3.822 2.496 4.135