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.
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(15859.4,0,15258.9,0,15498.6,0,15106.5,0,15023.6,0,12083,0,15761.3,0,16942.6,0,15070.3,0,13659.6,0,14768.9,0,14725.1,0,15998.1,0,15370.6,0,14956.9,0,15469.7,0,15101.8,0,11703.7,0,16283.6,0,16726.5,0,14968.9,0,14861,0,14583.3,0,15305.8,0,17903.9,0,16379.4,0,15420.3,0,17870.5,0,15912.8,0,13866.5,0,17823.2,0,17872,0,17422,0,16704.5,0,15991.2,0,16583.6,0,19123.5,0,17838.7,0,17209.4,0,18586.5,0,16258.1,0,15141.6,1,19202.1,1,17746.5,1,19090.1,1,18040.3,1,17515.5,1,17751.8,1,21072.4,1,17170,1,19439.5,1,19795.4,1,17574.9,1,16165.4,1,19464.6,1,19932.1,1,19961.2,1,17343.4,1,18924.2,1,18574.1,1,21350.6,1,18594.6,1,19823.1,1,20844.4,1,19640.2,1,17735.4,1,19813.6,1,22238.5,1,20682.2,1,17818.6,1,21872.1,1,22117,1,21865.9,1),dim=c(2,73),dimnames=list(c('uitvoer','dummy'),1:73))
> y <- array(NA,dim=c(2,73),dimnames=list(c('uitvoer','dummy'),1:73))
> 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 = '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
uitvoer dummy M1 M2 M3 M4 M5 M6 M7 M8 M9 M10 M11
1 15859.4 0 1 0 0 0 0 0 0 0 0 0 0
2 15258.9 0 0 1 0 0 0 0 0 0 0 0 0
3 15498.6 0 0 0 1 0 0 0 0 0 0 0 0
4 15106.5 0 0 0 0 1 0 0 0 0 0 0 0
5 15023.6 0 0 0 0 0 1 0 0 0 0 0 0
6 12083.0 0 0 0 0 0 0 1 0 0 0 0 0
7 15761.3 0 0 0 0 0 0 0 1 0 0 0 0
8 16942.6 0 0 0 0 0 0 0 0 1 0 0 0
9 15070.3 0 0 0 0 0 0 0 0 0 1 0 0
10 13659.6 0 0 0 0 0 0 0 0 0 0 1 0
11 14768.9 0 0 0 0 0 0 0 0 0 0 0 1
12 14725.1 0 0 0 0 0 0 0 0 0 0 0 0
13 15998.1 0 1 0 0 0 0 0 0 0 0 0 0
14 15370.6 0 0 1 0 0 0 0 0 0 0 0 0
15 14956.9 0 0 0 1 0 0 0 0 0 0 0 0
16 15469.7 0 0 0 0 1 0 0 0 0 0 0 0
17 15101.8 0 0 0 0 0 1 0 0 0 0 0 0
18 11703.7 0 0 0 0 0 0 1 0 0 0 0 0
19 16283.6 0 0 0 0 0 0 0 1 0 0 0 0
20 16726.5 0 0 0 0 0 0 0 0 1 0 0 0
21 14968.9 0 0 0 0 0 0 0 0 0 1 0 0
22 14861.0 0 0 0 0 0 0 0 0 0 0 1 0
23 14583.3 0 0 0 0 0 0 0 0 0 0 0 1
24 15305.8 0 0 0 0 0 0 0 0 0 0 0 0
25 17903.9 0 1 0 0 0 0 0 0 0 0 0 0
26 16379.4 0 0 1 0 0 0 0 0 0 0 0 0
27 15420.3 0 0 0 1 0 0 0 0 0 0 0 0
28 17870.5 0 0 0 0 1 0 0 0 0 0 0 0
29 15912.8 0 0 0 0 0 1 0 0 0 0 0 0
30 13866.5 0 0 0 0 0 0 1 0 0 0 0 0
31 17823.2 0 0 0 0 0 0 0 1 0 0 0 0
32 17872.0 0 0 0 0 0 0 0 0 1 0 0 0
33 17422.0 0 0 0 0 0 0 0 0 0 1 0 0
34 16704.5 0 0 0 0 0 0 0 0 0 0 1 0
35 15991.2 0 0 0 0 0 0 0 0 0 0 0 1
36 16583.6 0 0 0 0 0 0 0 0 0 0 0 0
37 19123.5 0 1 0 0 0 0 0 0 0 0 0 0
38 17838.7 0 0 1 0 0 0 0 0 0 0 0 0
39 17209.4 0 0 0 1 0 0 0 0 0 0 0 0
40 18586.5 0 0 0 0 1 0 0 0 0 0 0 0
41 16258.1 0 0 0 0 0 1 0 0 0 0 0 0
42 15141.6 1 0 0 0 0 0 1 0 0 0 0 0
43 19202.1 1 0 0 0 0 0 0 1 0 0 0 0
44 17746.5 1 0 0 0 0 0 0 0 1 0 0 0
45 19090.1 1 0 0 0 0 0 0 0 0 1 0 0
46 18040.3 1 0 0 0 0 0 0 0 0 0 1 0
47 17515.5 1 0 0 0 0 0 0 0 0 0 0 1
48 17751.8 1 0 0 0 0 0 0 0 0 0 0 0
49 21072.4 1 1 0 0 0 0 0 0 0 0 0 0
50 17170.0 1 0 1 0 0 0 0 0 0 0 0 0
51 19439.5 1 0 0 1 0 0 0 0 0 0 0 0
52 19795.4 1 0 0 0 1 0 0 0 0 0 0 0
53 17574.9 1 0 0 0 0 1 0 0 0 0 0 0
54 16165.4 1 0 0 0 0 0 1 0 0 0 0 0
55 19464.6 1 0 0 0 0 0 0 1 0 0 0 0
56 19932.1 1 0 0 0 0 0 0 0 1 0 0 0
57 19961.2 1 0 0 0 0 0 0 0 0 1 0 0
58 17343.4 1 0 0 0 0 0 0 0 0 0 1 0
59 18924.2 1 0 0 0 0 0 0 0 0 0 0 1
60 18574.1 1 0 0 0 0 0 0 0 0 0 0 0
61 21350.6 1 1 0 0 0 0 0 0 0 0 0 0
62 18594.6 1 0 1 0 0 0 0 0 0 0 0 0
63 19823.1 1 0 0 1 0 0 0 0 0 0 0 0
64 20844.4 1 0 0 0 1 0 0 0 0 0 0 0
65 19640.2 1 0 0 0 0 1 0 0 0 0 0 0
66 17735.4 1 0 0 0 0 0 1 0 0 0 0 0
67 19813.6 1 0 0 0 0 0 0 1 0 0 0 0
68 22238.5 1 0 0 0 0 0 0 0 1 0 0 0
69 20682.2 1 0 0 0 0 0 0 0 0 1 0 0
70 17818.6 1 0 0 0 0 0 0 0 0 0 1 0
71 21872.1 1 0 0 0 0 0 0 0 0 0 0 1
72 22117.0 1 0 0 0 0 0 0 0 0 0 0 0
73 21865.9 1 1 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) dummy M1 M2 M3 M4
15797.3 3424.6 1759.9 -170.1 119.2 1006.7
M5 M6 M7 M8 M9 M10
-353.6 -3060.3 548.5 1066.8 356.2 -1105.0
M11
-233.7
> (mysum <- summary(mylm))
Call:
lm(formula = df)
Residuals:
Min 1Q Median 3Q Max
-2542.2 -773.5 -137.6 771.9 2895.1
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 15797.3 526.7 29.994 < 2e-16 ***
dummy 3424.6 295.6 11.585 < 2e-16 ***
M1 1759.9 689.2 2.553 0.0132 *
M2 -170.1 716.6 -0.237 0.8132
M3 119.2 716.6 0.166 0.8685
M4 1006.7 716.6 1.405 0.1652
M5 -353.6 716.6 -0.493 0.6235
M6 -3060.3 714.9 -4.281 6.82e-05 ***
M7 548.5 714.9 0.767 0.4460
M8 1066.8 714.9 1.492 0.1409
M9 356.2 714.9 0.498 0.6201
M10 -1105.0 714.9 -1.546 0.1274
M11 -233.7 714.9 -0.327 0.7449
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 1238 on 60 degrees of freedom
Multiple R-squared: 0.7672, Adjusted R-squared: 0.7207
F-statistic: 16.48 on 12 and 60 DF, p-value: 9.258e-15
> 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.551188e-02 3.102377e-02 0.9844881
[2,] 2.912655e-03 5.825310e-03 0.9970873
[3,] 1.026960e-03 2.053921e-03 0.9989730
[4,] 5.336311e-04 1.067262e-03 0.9994664
[5,] 1.251732e-04 2.503464e-04 0.9998748
[6,] 3.214281e-05 6.428563e-05 0.9999679
[7,] 3.109296e-04 6.218592e-04 0.9996891
[8,] 1.213156e-04 2.426312e-04 0.9998787
[9,] 7.674175e-05 1.534835e-04 0.9999233
[10,] 5.158446e-03 1.031689e-02 0.9948416
[11,] 5.297542e-03 1.059508e-02 0.9947025
[12,] 3.354730e-03 6.709461e-03 0.9966453
[13,] 4.028160e-02 8.056321e-02 0.9597184
[14,] 3.029159e-02 6.058317e-02 0.9697084
[15,] 5.284764e-02 1.056953e-01 0.9471524
[16,] 6.768472e-02 1.353694e-01 0.9323153
[17,] 5.366675e-02 1.073335e-01 0.9463332
[18,] 9.348777e-02 1.869755e-01 0.9065122
[19,] 1.461121e-01 2.922242e-01 0.8538879
[20,] 1.399421e-01 2.798843e-01 0.8600579
[21,] 1.364609e-01 2.729218e-01 0.8635391
[22,] 1.879435e-01 3.758870e-01 0.8120565
[23,] 2.406882e-01 4.813763e-01 0.7593118
[24,] 2.344153e-01 4.688305e-01 0.7655847
[25,] 2.563655e-01 5.127310e-01 0.7436345
[26,] 2.030968e-01 4.061936e-01 0.7969032
[27,] 1.814725e-01 3.629450e-01 0.8185275
[28,] 1.321753e-01 2.643506e-01 0.8678247
[29,] 2.407508e-01 4.815016e-01 0.7592492
[30,] 2.143551e-01 4.287102e-01 0.7856449
[31,] 1.615429e-01 3.230859e-01 0.8384571
[32,] 2.305639e-01 4.611279e-01 0.7694361
[33,] 3.062072e-01 6.124143e-01 0.6937928
[34,] 2.559365e-01 5.118730e-01 0.7440635
[35,] 2.393803e-01 4.787605e-01 0.7606197
[36,] 1.824084e-01 3.648167e-01 0.8175916
[37,] 1.358254e-01 2.716508e-01 0.8641746
[38,] 1.339464e-01 2.678927e-01 0.8660536
[39,] 1.137675e-01 2.275349e-01 0.8862325
[40,] 6.606322e-02 1.321264e-01 0.9339368
[41,] 7.244288e-02 1.448858e-01 0.9275571
[42,] 4.031340e-02 8.062680e-02 0.9596866
> postscript(file="/var/www/html/rcomp/tmp/1mydf1230632192.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/2ds0q1230632192.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/36a261230632192.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/4supe1230632192.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/55cfi1230632192.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 = 73
Frequency = 1
1 2 3 4 5 6
-1697.743148 -368.266893 -417.833559 -1697.466893 -420.100226 -653.967006
7 8 9 10 11 12
-584.467006 78.532994 -1083.183673 -1032.667006 -794.667006 -1072.167006
13 14 15 16 17 18
-1559.043148 -256.566893 -959.533559 -1334.266893 -341.900226 -1033.267006
19 20 21 22 23 24
-62.167006 -137.567006 -1184.583673 168.732994 -980.267006 -491.467006
25 26 27 28 29 30
346.756852 752.233107 -496.133559 1066.533107 469.099774 1129.532994
31 32 33 34 35 36
1477.432994 1007.932994 1268.516327 2012.232994 427.632994 786.332994
37 38 39 40 41 42
1566.356852 2211.533107 1292.966441 1782.533107 814.399774 -1019.966327
43 44 45 46 47 48
-568.266327 -2542.166327 -487.982994 -76.566327 -1472.666327 -1470.066327
49 50 51 52 53 54
90.657531 -1881.766214 98.467119 -433.166214 -1293.399548 3.833673
55 56 57 58 59 60
-305.766327 -356.566327 383.117006 -773.466327 -63.966327 -647.766327
61 62 63 64 65 66
368.857531 -457.166214 482.067119 615.833786 771.900452 1573.833673
67 68 69 70 71 72
43.233673 1949.833673 1104.117006 -298.266327 2883.933673 2895.133673
73
884.157531
> postscript(file="/var/www/html/rcomp/tmp/6dal21230632192.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 = 73
Frequency = 1
lag(myerror, k = 1) myerror
0 -1697.743148 NA
1 -368.266893 -1697.743148
2 -417.833559 -368.266893
3 -1697.466893 -417.833559
4 -420.100226 -1697.466893
5 -653.967006 -420.100226
6 -584.467006 -653.967006
7 78.532994 -584.467006
8 -1083.183673 78.532994
9 -1032.667006 -1083.183673
10 -794.667006 -1032.667006
11 -1072.167006 -794.667006
12 -1559.043148 -1072.167006
13 -256.566893 -1559.043148
14 -959.533559 -256.566893
15 -1334.266893 -959.533559
16 -341.900226 -1334.266893
17 -1033.267006 -341.900226
18 -62.167006 -1033.267006
19 -137.567006 -62.167006
20 -1184.583673 -137.567006
21 168.732994 -1184.583673
22 -980.267006 168.732994
23 -491.467006 -980.267006
24 346.756852 -491.467006
25 752.233107 346.756852
26 -496.133559 752.233107
27 1066.533107 -496.133559
28 469.099774 1066.533107
29 1129.532994 469.099774
30 1477.432994 1129.532994
31 1007.932994 1477.432994
32 1268.516327 1007.932994
33 2012.232994 1268.516327
34 427.632994 2012.232994
35 786.332994 427.632994
36 1566.356852 786.332994
37 2211.533107 1566.356852
38 1292.966441 2211.533107
39 1782.533107 1292.966441
40 814.399774 1782.533107
41 -1019.966327 814.399774
42 -568.266327 -1019.966327
43 -2542.166327 -568.266327
44 -487.982994 -2542.166327
45 -76.566327 -487.982994
46 -1472.666327 -76.566327
47 -1470.066327 -1472.666327
48 90.657531 -1470.066327
49 -1881.766214 90.657531
50 98.467119 -1881.766214
51 -433.166214 98.467119
52 -1293.399548 -433.166214
53 3.833673 -1293.399548
54 -305.766327 3.833673
55 -356.566327 -305.766327
56 383.117006 -356.566327
57 -773.466327 383.117006
58 -63.966327 -773.466327
59 -647.766327 -63.966327
60 368.857531 -647.766327
61 -457.166214 368.857531
62 482.067119 -457.166214
63 615.833786 482.067119
64 771.900452 615.833786
65 1573.833673 771.900452
66 43.233673 1573.833673
67 1949.833673 43.233673
68 1104.117006 1949.833673
69 -298.266327 1104.117006
70 2883.933673 -298.266327
71 2895.133673 2883.933673
72 884.157531 2895.133673
73 NA 884.157531
> dum1 <- dum[2:length(myerror),]
> dum1
lag(myerror, k = 1) myerror
[1,] -368.266893 -1697.743148
[2,] -417.833559 -368.266893
[3,] -1697.466893 -417.833559
[4,] -420.100226 -1697.466893
[5,] -653.967006 -420.100226
[6,] -584.467006 -653.967006
[7,] 78.532994 -584.467006
[8,] -1083.183673 78.532994
[9,] -1032.667006 -1083.183673
[10,] -794.667006 -1032.667006
[11,] -1072.167006 -794.667006
[12,] -1559.043148 -1072.167006
[13,] -256.566893 -1559.043148
[14,] -959.533559 -256.566893
[15,] -1334.266893 -959.533559
[16,] -341.900226 -1334.266893
[17,] -1033.267006 -341.900226
[18,] -62.167006 -1033.267006
[19,] -137.567006 -62.167006
[20,] -1184.583673 -137.567006
[21,] 168.732994 -1184.583673
[22,] -980.267006 168.732994
[23,] -491.467006 -980.267006
[24,] 346.756852 -491.467006
[25,] 752.233107 346.756852
[26,] -496.133559 752.233107
[27,] 1066.533107 -496.133559
[28,] 469.099774 1066.533107
[29,] 1129.532994 469.099774
[30,] 1477.432994 1129.532994
[31,] 1007.932994 1477.432994
[32,] 1268.516327 1007.932994
[33,] 2012.232994 1268.516327
[34,] 427.632994 2012.232994
[35,] 786.332994 427.632994
[36,] 1566.356852 786.332994
[37,] 2211.533107 1566.356852
[38,] 1292.966441 2211.533107
[39,] 1782.533107 1292.966441
[40,] 814.399774 1782.533107
[41,] -1019.966327 814.399774
[42,] -568.266327 -1019.966327
[43,] -2542.166327 -568.266327
[44,] -487.982994 -2542.166327
[45,] -76.566327 -487.982994
[46,] -1472.666327 -76.566327
[47,] -1470.066327 -1472.666327
[48,] 90.657531 -1470.066327
[49,] -1881.766214 90.657531
[50,] 98.467119 -1881.766214
[51,] -433.166214 98.467119
[52,] -1293.399548 -433.166214
[53,] 3.833673 -1293.399548
[54,] -305.766327 3.833673
[55,] -356.566327 -305.766327
[56,] 383.117006 -356.566327
[57,] -773.466327 383.117006
[58,] -63.966327 -773.466327
[59,] -647.766327 -63.966327
[60,] 368.857531 -647.766327
[61,] -457.166214 368.857531
[62,] 482.067119 -457.166214
[63,] 615.833786 482.067119
[64,] 771.900452 615.833786
[65,] 1573.833673 771.900452
[66,] 43.233673 1573.833673
[67,] 1949.833673 43.233673
[68,] 1104.117006 1949.833673
[69,] -298.266327 1104.117006
[70,] 2883.933673 -298.266327
[71,] 2895.133673 2883.933673
[72,] 884.157531 2895.133673
> z <- as.data.frame(dum1)
> z
lag(myerror, k = 1) myerror
1 -368.266893 -1697.743148
2 -417.833559 -368.266893
3 -1697.466893 -417.833559
4 -420.100226 -1697.466893
5 -653.967006 -420.100226
6 -584.467006 -653.967006
7 78.532994 -584.467006
8 -1083.183673 78.532994
9 -1032.667006 -1083.183673
10 -794.667006 -1032.667006
11 -1072.167006 -794.667006
12 -1559.043148 -1072.167006
13 -256.566893 -1559.043148
14 -959.533559 -256.566893
15 -1334.266893 -959.533559
16 -341.900226 -1334.266893
17 -1033.267006 -341.900226
18 -62.167006 -1033.267006
19 -137.567006 -62.167006
20 -1184.583673 -137.567006
21 168.732994 -1184.583673
22 -980.267006 168.732994
23 -491.467006 -980.267006
24 346.756852 -491.467006
25 752.233107 346.756852
26 -496.133559 752.233107
27 1066.533107 -496.133559
28 469.099774 1066.533107
29 1129.532994 469.099774
30 1477.432994 1129.532994
31 1007.932994 1477.432994
32 1268.516327 1007.932994
33 2012.232994 1268.516327
34 427.632994 2012.232994
35 786.332994 427.632994
36 1566.356852 786.332994
37 2211.533107 1566.356852
38 1292.966441 2211.533107
39 1782.533107 1292.966441
40 814.399774 1782.533107
41 -1019.966327 814.399774
42 -568.266327 -1019.966327
43 -2542.166327 -568.266327
44 -487.982994 -2542.166327
45 -76.566327 -487.982994
46 -1472.666327 -76.566327
47 -1470.066327 -1472.666327
48 90.657531 -1470.066327
49 -1881.766214 90.657531
50 98.467119 -1881.766214
51 -433.166214 98.467119
52 -1293.399548 -433.166214
53 3.833673 -1293.399548
54 -305.766327 3.833673
55 -356.566327 -305.766327
56 383.117006 -356.566327
57 -773.466327 383.117006
58 -63.966327 -773.466327
59 -647.766327 -63.966327
60 368.857531 -647.766327
61 -457.166214 368.857531
62 482.067119 -457.166214
63 615.833786 482.067119
64 771.900452 615.833786
65 1573.833673 771.900452
66 43.233673 1573.833673
67 1949.833673 43.233673
68 1104.117006 1949.833673
69 -298.266327 1104.117006
70 2883.933673 -298.266327
71 2895.133673 2883.933673
72 884.157531 2895.133673
> 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/7426v1230632192.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/8q0et1230632192.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/9859c1230632192.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/10j0k71230632192.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/11fynu1230632192.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/129rn71230632192.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/13qsnp1230632193.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/14hs7d1230632193.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/159iwe1230632193.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/16i4jc1230632193.tab")
+ }
>
> system("convert tmp/1mydf1230632192.ps tmp/1mydf1230632192.png")
> system("convert tmp/2ds0q1230632192.ps tmp/2ds0q1230632192.png")
> system("convert tmp/36a261230632192.ps tmp/36a261230632192.png")
> system("convert tmp/4supe1230632192.ps tmp/4supe1230632192.png")
> system("convert tmp/55cfi1230632192.ps tmp/55cfi1230632192.png")
> system("convert tmp/6dal21230632192.ps tmp/6dal21230632192.png")
> system("convert tmp/7426v1230632192.ps tmp/7426v1230632192.png")
> system("convert tmp/8q0et1230632192.ps tmp/8q0et1230632192.png")
> system("convert tmp/9859c1230632192.ps tmp/9859c1230632192.png")
> system("convert tmp/10j0k71230632192.ps tmp/10j0k71230632192.png")
>
>
> proc.time()
user system elapsed
2.718 1.624 3.516