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(46.85,48.05,54.63,53.22,49.87,56.42,59.03,64.99,65.55,62.27,58.34,59.45,65.54,61.93,62.97,70.16,70.96,70.97,74.46,73.08,63.90,59.14,59.40,62.09,54.35,59.39,60.74,64.04,63.53,67.53,74.15,72.36,79.63,85.66,94.63,91.74,92.93,95.35,105.42,112.46,125.46,134.02,133.48,116.69,103.76,76.72,57.44,42.04,41.92,39.26,48.06,49.95,59.21,69.70,64.29,71.14,69.47,75.82,78.15,74.60),dim=c(1,60),dimnames=list(c('Crudeoilprice'),1:60))
> y <- array(NA,dim=c(1,60),dimnames=list(c('Crudeoilprice'),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
Crudeoilprice M1 M2 M3 M4 M5 M6 M7 M8 M9 M10 M11 t
1 46.85 1 0 0 0 0 0 0 0 0 0 0 1
2 48.05 0 1 0 0 0 0 0 0 0 0 0 2
3 54.63 0 0 1 0 0 0 0 0 0 0 0 3
4 53.22 0 0 0 1 0 0 0 0 0 0 0 4
5 49.87 0 0 0 0 1 0 0 0 0 0 0 5
6 56.42 0 0 0 0 0 1 0 0 0 0 0 6
7 59.03 0 0 0 0 0 0 1 0 0 0 0 7
8 64.99 0 0 0 0 0 0 0 1 0 0 0 8
9 65.55 0 0 0 0 0 0 0 0 1 0 0 9
10 62.27 0 0 0 0 0 0 0 0 0 1 0 10
11 58.34 0 0 0 0 0 0 0 0 0 0 1 11
12 59.45 0 0 0 0 0 0 0 0 0 0 0 12
13 65.54 1 0 0 0 0 0 0 0 0 0 0 13
14 61.93 0 1 0 0 0 0 0 0 0 0 0 14
15 62.97 0 0 1 0 0 0 0 0 0 0 0 15
16 70.16 0 0 0 1 0 0 0 0 0 0 0 16
17 70.96 0 0 0 0 1 0 0 0 0 0 0 17
18 70.97 0 0 0 0 0 1 0 0 0 0 0 18
19 74.46 0 0 0 0 0 0 1 0 0 0 0 19
20 73.08 0 0 0 0 0 0 0 1 0 0 0 20
21 63.90 0 0 0 0 0 0 0 0 1 0 0 21
22 59.14 0 0 0 0 0 0 0 0 0 1 0 22
23 59.40 0 0 0 0 0 0 0 0 0 0 1 23
24 62.09 0 0 0 0 0 0 0 0 0 0 0 24
25 54.35 1 0 0 0 0 0 0 0 0 0 0 25
26 59.39 0 1 0 0 0 0 0 0 0 0 0 26
27 60.74 0 0 1 0 0 0 0 0 0 0 0 27
28 64.04 0 0 0 1 0 0 0 0 0 0 0 28
29 63.53 0 0 0 0 1 0 0 0 0 0 0 29
30 67.53 0 0 0 0 0 1 0 0 0 0 0 30
31 74.15 0 0 0 0 0 0 1 0 0 0 0 31
32 72.36 0 0 0 0 0 0 0 1 0 0 0 32
33 79.63 0 0 0 0 0 0 0 0 1 0 0 33
34 85.66 0 0 0 0 0 0 0 0 0 1 0 34
35 94.63 0 0 0 0 0 0 0 0 0 0 1 35
36 91.74 0 0 0 0 0 0 0 0 0 0 0 36
37 92.93 1 0 0 0 0 0 0 0 0 0 0 37
38 95.35 0 1 0 0 0 0 0 0 0 0 0 38
39 105.42 0 0 1 0 0 0 0 0 0 0 0 39
40 112.46 0 0 0 1 0 0 0 0 0 0 0 40
41 125.46 0 0 0 0 1 0 0 0 0 0 0 41
42 134.02 0 0 0 0 0 1 0 0 0 0 0 42
43 133.48 0 0 0 0 0 0 1 0 0 0 0 43
44 116.69 0 0 0 0 0 0 0 1 0 0 0 44
45 103.76 0 0 0 0 0 0 0 0 1 0 0 45
46 76.72 0 0 0 0 0 0 0 0 0 1 0 46
47 57.44 0 0 0 0 0 0 0 0 0 0 1 47
48 42.04 0 0 0 0 0 0 0 0 0 0 0 48
49 41.92 1 0 0 0 0 0 0 0 0 0 0 49
50 39.26 0 1 0 0 0 0 0 0 0 0 0 50
51 48.06 0 0 1 0 0 0 0 0 0 0 0 51
52 49.95 0 0 0 1 0 0 0 0 0 0 0 52
53 59.21 0 0 0 0 1 0 0 0 0 0 0 53
54 69.70 0 0 0 0 0 1 0 0 0 0 0 54
55 64.29 0 0 0 0 0 0 1 0 0 0 0 55
56 71.14 0 0 0 0 0 0 0 1 0 0 0 56
57 69.47 0 0 0 0 0 0 0 0 1 0 0 57
58 75.82 0 0 0 0 0 0 0 0 0 1 0 58
59 78.15 0 0 0 0 0 0 0 0 0 0 1 59
60 74.60 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) M1 M2 M3 M4 M5
52.8098 -1.6405 -1.5285 3.6736 6.9096 10.3837
M6 M7 M8 M9 M10 M11
15.9397 16.9278 15.1318 11.5759 6.6699 3.9740
t
0.3660
> (mysum <- summary(mylm))
Call:
lm(formula = df)
Residuals:
Min 1Q Median 3Q Max
-30.319 -10.757 -3.457 3.522 49.901
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 52.8098 11.4678 4.605 3.15e-05 ***
M1 -1.6405 13.9513 -0.118 0.9069
M2 -1.5285 13.9304 -0.110 0.9131
M3 3.6736 13.9115 0.264 0.7929
M4 6.9096 13.8946 0.497 0.6213
M5 10.3837 13.8797 0.748 0.4581
M6 15.9397 13.8667 1.149 0.2562
M7 16.9278 13.8557 1.222 0.2279
M8 15.1318 13.8467 1.093 0.2800
M9 11.5759 13.8397 0.836 0.4072
M10 6.6699 13.8347 0.482 0.6320
M11 3.9740 13.8317 0.287 0.7751
t 0.3660 0.1664 2.200 0.0328 *
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 21.87 on 47 degrees of freedom
Multiple R-squared: 0.1851, Adjusted R-squared: -0.023
F-statistic: 0.8895 on 12 and 47 DF, p-value: 0.5633
> 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.630539e-03 5.261078e-03 0.9973695
[2,] 6.318686e-04 1.263737e-03 0.9993681
[3,] 7.346417e-05 1.469283e-04 0.9999265
[4,] 7.567463e-06 1.513493e-05 0.9999924
[5,] 2.889956e-06 5.779912e-06 0.9999971
[6,] 1.374056e-05 2.748111e-05 0.9999863
[7,] 1.822256e-05 3.644511e-05 0.9999818
[8,] 7.525752e-06 1.505150e-05 0.9999925
[9,] 2.227260e-06 4.454519e-06 0.9999978
[10,] 2.761033e-06 5.522067e-06 0.9999972
[11,] 8.206033e-07 1.641207e-06 0.9999992
[12,] 3.061504e-07 6.123009e-07 0.9999997
[13,] 1.051969e-07 2.103939e-07 0.9999999
[14,] 4.607770e-08 9.215540e-08 1.0000000
[15,] 3.041727e-08 6.083454e-08 1.0000000
[16,] 1.921526e-08 3.843053e-08 1.0000000
[17,] 2.113928e-08 4.227857e-08 1.0000000
[18,] 3.030144e-08 6.060288e-08 1.0000000
[19,] 1.606630e-07 3.213261e-07 0.9999998
[20,] 2.625107e-06 5.250215e-06 0.9999974
[21,] 5.255956e-06 1.051191e-05 0.9999947
[22,] 8.525530e-06 1.705106e-05 0.9999915
[23,] 1.188807e-05 2.377614e-05 0.9999881
[24,] 2.722371e-05 5.444742e-05 0.9999728
[25,] 7.766537e-05 1.553307e-04 0.9999223
[26,] 7.142055e-04 1.428411e-03 0.9992858
[27,] 4.879081e-03 9.758161e-03 0.9951209
[28,] 4.443978e-02 8.887956e-02 0.9555602
[29,] 1.183999e-01 2.367997e-01 0.8816001
> postscript(file="/var/www/html/rcomp/tmp/1e1qm1291149381.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/2osp71291149381.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/3osp71291149381.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/4osp71291149381.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/5z2ps1291149381.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 = 60
Frequency = 1
1 2 3 4 5 6
-4.6851667 -3.9631667 -2.9511667 -7.9631667 -15.1531667 -14.5251667
7 8 9 10 11 12
-13.2691667 -5.8791667 -2.1291667 -0.8691667 -2.4691667 2.2488333
13 14 15 16 17 18
9.6134167 5.5254167 0.9974167 4.5854167 1.5454167 -4.3665833
19 20 21 22 23 24
-2.2305833 -2.1805833 -8.1705833 -8.3905833 -5.8005833 0.4974167
25 26 27 28 29 30
-5.9680000 -1.4060000 -5.6240000 -5.9260000 -10.2760000 -12.1980000
31 32 33 34 35 36
-6.9320000 -7.2920000 3.1680000 13.7380000 25.0380000 25.7560000
37 38 39 40 41 42
28.2205833 30.1625833 34.6645833 38.1025833 47.2625833 49.9005833
43 44 45 46 47 48
48.0065833 32.6465833 22.9065833 0.4065833 -16.5434167 -28.3354167
49 50 51 52 53 54
-27.1808333 -30.3188333 -27.0868333 -28.7988333 -23.3788333 -18.8108333
55 56 57 58 59 60
-25.5748333 -17.2948333 -15.7748333 -4.8848333 -0.2248333 -0.1668333
> postscript(file="/var/www/html/rcomp/tmp/6z2ps1291149381.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 = 60
Frequency = 1
lag(myerror, k = 1) myerror
0 -4.6851667 NA
1 -3.9631667 -4.6851667
2 -2.9511667 -3.9631667
3 -7.9631667 -2.9511667
4 -15.1531667 -7.9631667
5 -14.5251667 -15.1531667
6 -13.2691667 -14.5251667
7 -5.8791667 -13.2691667
8 -2.1291667 -5.8791667
9 -0.8691667 -2.1291667
10 -2.4691667 -0.8691667
11 2.2488333 -2.4691667
12 9.6134167 2.2488333
13 5.5254167 9.6134167
14 0.9974167 5.5254167
15 4.5854167 0.9974167
16 1.5454167 4.5854167
17 -4.3665833 1.5454167
18 -2.2305833 -4.3665833
19 -2.1805833 -2.2305833
20 -8.1705833 -2.1805833
21 -8.3905833 -8.1705833
22 -5.8005833 -8.3905833
23 0.4974167 -5.8005833
24 -5.9680000 0.4974167
25 -1.4060000 -5.9680000
26 -5.6240000 -1.4060000
27 -5.9260000 -5.6240000
28 -10.2760000 -5.9260000
29 -12.1980000 -10.2760000
30 -6.9320000 -12.1980000
31 -7.2920000 -6.9320000
32 3.1680000 -7.2920000
33 13.7380000 3.1680000
34 25.0380000 13.7380000
35 25.7560000 25.0380000
36 28.2205833 25.7560000
37 30.1625833 28.2205833
38 34.6645833 30.1625833
39 38.1025833 34.6645833
40 47.2625833 38.1025833
41 49.9005833 47.2625833
42 48.0065833 49.9005833
43 32.6465833 48.0065833
44 22.9065833 32.6465833
45 0.4065833 22.9065833
46 -16.5434167 0.4065833
47 -28.3354167 -16.5434167
48 -27.1808333 -28.3354167
49 -30.3188333 -27.1808333
50 -27.0868333 -30.3188333
51 -28.7988333 -27.0868333
52 -23.3788333 -28.7988333
53 -18.8108333 -23.3788333
54 -25.5748333 -18.8108333
55 -17.2948333 -25.5748333
56 -15.7748333 -17.2948333
57 -4.8848333 -15.7748333
58 -0.2248333 -4.8848333
59 -0.1668333 -0.2248333
60 NA -0.1668333
> dum1 <- dum[2:length(myerror),]
> dum1
lag(myerror, k = 1) myerror
[1,] -3.9631667 -4.6851667
[2,] -2.9511667 -3.9631667
[3,] -7.9631667 -2.9511667
[4,] -15.1531667 -7.9631667
[5,] -14.5251667 -15.1531667
[6,] -13.2691667 -14.5251667
[7,] -5.8791667 -13.2691667
[8,] -2.1291667 -5.8791667
[9,] -0.8691667 -2.1291667
[10,] -2.4691667 -0.8691667
[11,] 2.2488333 -2.4691667
[12,] 9.6134167 2.2488333
[13,] 5.5254167 9.6134167
[14,] 0.9974167 5.5254167
[15,] 4.5854167 0.9974167
[16,] 1.5454167 4.5854167
[17,] -4.3665833 1.5454167
[18,] -2.2305833 -4.3665833
[19,] -2.1805833 -2.2305833
[20,] -8.1705833 -2.1805833
[21,] -8.3905833 -8.1705833
[22,] -5.8005833 -8.3905833
[23,] 0.4974167 -5.8005833
[24,] -5.9680000 0.4974167
[25,] -1.4060000 -5.9680000
[26,] -5.6240000 -1.4060000
[27,] -5.9260000 -5.6240000
[28,] -10.2760000 -5.9260000
[29,] -12.1980000 -10.2760000
[30,] -6.9320000 -12.1980000
[31,] -7.2920000 -6.9320000
[32,] 3.1680000 -7.2920000
[33,] 13.7380000 3.1680000
[34,] 25.0380000 13.7380000
[35,] 25.7560000 25.0380000
[36,] 28.2205833 25.7560000
[37,] 30.1625833 28.2205833
[38,] 34.6645833 30.1625833
[39,] 38.1025833 34.6645833
[40,] 47.2625833 38.1025833
[41,] 49.9005833 47.2625833
[42,] 48.0065833 49.9005833
[43,] 32.6465833 48.0065833
[44,] 22.9065833 32.6465833
[45,] 0.4065833 22.9065833
[46,] -16.5434167 0.4065833
[47,] -28.3354167 -16.5434167
[48,] -27.1808333 -28.3354167
[49,] -30.3188333 -27.1808333
[50,] -27.0868333 -30.3188333
[51,] -28.7988333 -27.0868333
[52,] -23.3788333 -28.7988333
[53,] -18.8108333 -23.3788333
[54,] -25.5748333 -18.8108333
[55,] -17.2948333 -25.5748333
[56,] -15.7748333 -17.2948333
[57,] -4.8848333 -15.7748333
[58,] -0.2248333 -4.8848333
[59,] -0.1668333 -0.2248333
> z <- as.data.frame(dum1)
> z
lag(myerror, k = 1) myerror
1 -3.9631667 -4.6851667
2 -2.9511667 -3.9631667
3 -7.9631667 -2.9511667
4 -15.1531667 -7.9631667
5 -14.5251667 -15.1531667
6 -13.2691667 -14.5251667
7 -5.8791667 -13.2691667
8 -2.1291667 -5.8791667
9 -0.8691667 -2.1291667
10 -2.4691667 -0.8691667
11 2.2488333 -2.4691667
12 9.6134167 2.2488333
13 5.5254167 9.6134167
14 0.9974167 5.5254167
15 4.5854167 0.9974167
16 1.5454167 4.5854167
17 -4.3665833 1.5454167
18 -2.2305833 -4.3665833
19 -2.1805833 -2.2305833
20 -8.1705833 -2.1805833
21 -8.3905833 -8.1705833
22 -5.8005833 -8.3905833
23 0.4974167 -5.8005833
24 -5.9680000 0.4974167
25 -1.4060000 -5.9680000
26 -5.6240000 -1.4060000
27 -5.9260000 -5.6240000
28 -10.2760000 -5.9260000
29 -12.1980000 -10.2760000
30 -6.9320000 -12.1980000
31 -7.2920000 -6.9320000
32 3.1680000 -7.2920000
33 13.7380000 3.1680000
34 25.0380000 13.7380000
35 25.7560000 25.0380000
36 28.2205833 25.7560000
37 30.1625833 28.2205833
38 34.6645833 30.1625833
39 38.1025833 34.6645833
40 47.2625833 38.1025833
41 49.9005833 47.2625833
42 48.0065833 49.9005833
43 32.6465833 48.0065833
44 22.9065833 32.6465833
45 0.4065833 22.9065833
46 -16.5434167 0.4065833
47 -28.3354167 -16.5434167
48 -27.1808333 -28.3354167
49 -30.3188333 -27.1808333
50 -27.0868333 -30.3188333
51 -28.7988333 -27.0868333
52 -23.3788333 -28.7988333
53 -18.8108333 -23.3788333
54 -25.5748333 -18.8108333
55 -17.2948333 -25.5748333
56 -15.7748333 -17.2948333
57 -4.8848333 -15.7748333
58 -0.2248333 -4.8848333
59 -0.1668333 -0.2248333
> 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/7sb6d1291149381.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/8sb6d1291149381.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/922ng1291149381.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/1022ng1291149381.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/11o3441291149381.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/1293291291149381.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/13y4z31291149381.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/148ego1291149381.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/15cexu1291149381.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/16xwd01291149381.tab")
+ }
>
> try(system("convert tmp/1e1qm1291149381.ps tmp/1e1qm1291149381.png",intern=TRUE))
character(0)
> try(system("convert tmp/2osp71291149381.ps tmp/2osp71291149381.png",intern=TRUE))
character(0)
> try(system("convert tmp/3osp71291149381.ps tmp/3osp71291149381.png",intern=TRUE))
character(0)
> try(system("convert tmp/4osp71291149381.ps tmp/4osp71291149381.png",intern=TRUE))
character(0)
> try(system("convert tmp/5z2ps1291149381.ps tmp/5z2ps1291149381.png",intern=TRUE))
character(0)
> try(system("convert tmp/6z2ps1291149381.ps tmp/6z2ps1291149381.png",intern=TRUE))
character(0)
> try(system("convert tmp/7sb6d1291149381.ps tmp/7sb6d1291149381.png",intern=TRUE))
character(0)
> try(system("convert tmp/8sb6d1291149381.ps tmp/8sb6d1291149381.png",intern=TRUE))
character(0)
> try(system("convert tmp/922ng1291149381.ps tmp/922ng1291149381.png",intern=TRUE))
character(0)
> try(system("convert tmp/1022ng1291149381.ps tmp/1022ng1291149381.png",intern=TRUE))
character(0)
>
>
> proc.time()
user system elapsed
2.370 1.564 5.561