R version 2.7.0 (2008-04-22)
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(97.7,0,101.5,0,119.6,0,108.1,0,117.8,0,125.5,0,89.2,0,92.3,0,104.6,0,122.8,0,96.0,0,94.6,0,93.3,0,101.1,0,114.2,0,104.7,0,113.3,0,118.2,0,83.6,0,73.9,0,99.5,0,97.7,0,103.0,0,106.3,0,92.2,0,101.8,0,122.8,0,111.8,0,106.3,0,121.5,0,81.9,0,85.4,0,110.9,0,117.3,0,106.3,0,105.5,0,101.3,0,105.9,0,126.3,0,111.9,0,108.9,0,127.2,0,94.2,0,85.7,0,116.2,0,107.2,0,110.6,0,112.0,0,104.5,0,112.0,0,132.8,0,110.8,0,128.7,0,136.8,0,94.9,0,88.8,0,123.2,0,125.3,0,122.7,0,125.7,0,116.3,0,118.7,0,142.0,0,127.9,0,131.9,0,152.3,0,110.8,1,99.1,1,135.0,1,133.2,1,131.0,1,133.9,1,119.9,1,136.9,1,148.9,1,145.1,1,142.4,1,159.6,1,120.7,1,109.0,1,142.0,1),dim=c(2,81),dimnames=list(c('Y','X'),1:81))
> y <- array(NA,dim=c(2,81),dimnames=list(c('Y','X'),1:81))
> 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
Y X M1 M2 M3 M4 M5 M6 M7 M8 M9 M10 M11 t
1 97.7 0 1 0 0 0 0 0 0 0 0 0 0 1
2 101.5 0 0 1 0 0 0 0 0 0 0 0 0 2
3 119.6 0 0 0 1 0 0 0 0 0 0 0 0 3
4 108.1 0 0 0 0 1 0 0 0 0 0 0 0 4
5 117.8 0 0 0 0 0 1 0 0 0 0 0 0 5
6 125.5 0 0 0 0 0 0 1 0 0 0 0 0 6
7 89.2 0 0 0 0 0 0 0 1 0 0 0 0 7
8 92.3 0 0 0 0 0 0 0 0 1 0 0 0 8
9 104.6 0 0 0 0 0 0 0 0 0 1 0 0 9
10 122.8 0 0 0 0 0 0 0 0 0 0 1 0 10
11 96.0 0 0 0 0 0 0 0 0 0 0 0 1 11
12 94.6 0 0 0 0 0 0 0 0 0 0 0 0 12
13 93.3 0 1 0 0 0 0 0 0 0 0 0 0 13
14 101.1 0 0 1 0 0 0 0 0 0 0 0 0 14
15 114.2 0 0 0 1 0 0 0 0 0 0 0 0 15
16 104.7 0 0 0 0 1 0 0 0 0 0 0 0 16
17 113.3 0 0 0 0 0 1 0 0 0 0 0 0 17
18 118.2 0 0 0 0 0 0 1 0 0 0 0 0 18
19 83.6 0 0 0 0 0 0 0 1 0 0 0 0 19
20 73.9 0 0 0 0 0 0 0 0 1 0 0 0 20
21 99.5 0 0 0 0 0 0 0 0 0 1 0 0 21
22 97.7 0 0 0 0 0 0 0 0 0 0 1 0 22
23 103.0 0 0 0 0 0 0 0 0 0 0 0 1 23
24 106.3 0 0 0 0 0 0 0 0 0 0 0 0 24
25 92.2 0 1 0 0 0 0 0 0 0 0 0 0 25
26 101.8 0 0 1 0 0 0 0 0 0 0 0 0 26
27 122.8 0 0 0 1 0 0 0 0 0 0 0 0 27
28 111.8 0 0 0 0 1 0 0 0 0 0 0 0 28
29 106.3 0 0 0 0 0 1 0 0 0 0 0 0 29
30 121.5 0 0 0 0 0 0 1 0 0 0 0 0 30
31 81.9 0 0 0 0 0 0 0 1 0 0 0 0 31
32 85.4 0 0 0 0 0 0 0 0 1 0 0 0 32
33 110.9 0 0 0 0 0 0 0 0 0 1 0 0 33
34 117.3 0 0 0 0 0 0 0 0 0 0 1 0 34
35 106.3 0 0 0 0 0 0 0 0 0 0 0 1 35
36 105.5 0 0 0 0 0 0 0 0 0 0 0 0 36
37 101.3 0 1 0 0 0 0 0 0 0 0 0 0 37
38 105.9 0 0 1 0 0 0 0 0 0 0 0 0 38
39 126.3 0 0 0 1 0 0 0 0 0 0 0 0 39
40 111.9 0 0 0 0 1 0 0 0 0 0 0 0 40
41 108.9 0 0 0 0 0 1 0 0 0 0 0 0 41
42 127.2 0 0 0 0 0 0 1 0 0 0 0 0 42
43 94.2 0 0 0 0 0 0 0 1 0 0 0 0 43
44 85.7 0 0 0 0 0 0 0 0 1 0 0 0 44
45 116.2 0 0 0 0 0 0 0 0 0 1 0 0 45
46 107.2 0 0 0 0 0 0 0 0 0 0 1 0 46
47 110.6 0 0 0 0 0 0 0 0 0 0 0 1 47
48 112.0 0 0 0 0 0 0 0 0 0 0 0 0 48
49 104.5 0 1 0 0 0 0 0 0 0 0 0 0 49
50 112.0 0 0 1 0 0 0 0 0 0 0 0 0 50
51 132.8 0 0 0 1 0 0 0 0 0 0 0 0 51
52 110.8 0 0 0 0 1 0 0 0 0 0 0 0 52
53 128.7 0 0 0 0 0 1 0 0 0 0 0 0 53
54 136.8 0 0 0 0 0 0 1 0 0 0 0 0 54
55 94.9 0 0 0 0 0 0 0 1 0 0 0 0 55
56 88.8 0 0 0 0 0 0 0 0 1 0 0 0 56
57 123.2 0 0 0 0 0 0 0 0 0 1 0 0 57
58 125.3 0 0 0 0 0 0 0 0 0 0 1 0 58
59 122.7 0 0 0 0 0 0 0 0 0 0 0 1 59
60 125.7 0 0 0 0 0 0 0 0 0 0 0 0 60
61 116.3 0 1 0 0 0 0 0 0 0 0 0 0 61
62 118.7 0 0 1 0 0 0 0 0 0 0 0 0 62
63 142.0 0 0 0 1 0 0 0 0 0 0 0 0 63
64 127.9 0 0 0 0 1 0 0 0 0 0 0 0 64
65 131.9 0 0 0 0 0 1 0 0 0 0 0 0 65
66 152.3 0 0 0 0 0 0 1 0 0 0 0 0 66
67 110.8 1 0 0 0 0 0 0 1 0 0 0 0 67
68 99.1 1 0 0 0 0 0 0 0 1 0 0 0 68
69 135.0 1 0 0 0 0 0 0 0 0 1 0 0 69
70 133.2 1 0 0 0 0 0 0 0 0 0 1 0 70
71 131.0 1 0 0 0 0 0 0 0 0 0 0 1 71
72 133.9 1 0 0 0 0 0 0 0 0 0 0 0 72
73 119.9 1 1 0 0 0 0 0 0 0 0 0 0 73
74 136.9 1 0 1 0 0 0 0 0 0 0 0 0 74
75 148.9 1 0 0 1 0 0 0 0 0 0 0 0 75
76 145.1 1 0 0 0 1 0 0 0 0 0 0 0 76
77 142.4 1 0 0 0 0 1 0 0 0 0 0 0 77
78 159.6 1 0 0 0 0 0 1 0 0 0 0 0 78
79 120.7 1 0 0 0 0 0 0 1 0 0 0 0 79
80 109.0 1 0 0 0 0 0 0 0 1 0 0 0 80
81 142.0 1 0 0 0 0 0 0 0 0 1 0 0 81
> k <- length(x[1,])
> df <- as.data.frame(x)
> (mylm <- lm(df))
Call:
lm(formula = df)
Coefficients:
(Intercept) X M1 M2 M3 M4
97.0556 11.2265 -7.4573 -0.2638 17.7868 5.1232
M5 M6 M7 M8 M9 M10
8.9309 21.7102 -18.2001 -24.4066 3.4297 4.9202
M11 t
-1.0649 0.3351
> (mysum <- summary(mylm))
Call:
lm(formula = df)
Residuals:
Min 1Q Median 3Q Max
-11.6475 -3.6184 -0.6429 3.8897 17.4734
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 97.05562 2.90248 33.439 < 2e-16 ***
X 11.22651 2.41638 4.646 1.63e-05 ***
M1 -7.45731 3.44247 -2.166 0.0339 *
M2 -0.26382 3.44084 -0.077 0.9391
M3 17.78682 3.43967 5.171 2.29e-06 ***
M4 5.12317 3.43896 1.490 0.1410
M5 8.93095 3.43872 2.597 0.0115 *
M6 21.71015 3.43894 6.313 2.51e-08 ***
M7 -18.20014 3.44846 -5.278 1.52e-06 ***
M8 -24.40665 3.44692 -7.081 1.09e-09 ***
M9 3.42970 3.44584 0.995 0.3232
M10 4.92016 3.56914 1.379 0.1726
M11 -1.06492 3.56848 -0.298 0.7663
t 0.33508 0.03986 8.407 4.44e-12 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 6.18 on 67 degrees of freedom
Multiple R-squared: 0.8966, Adjusted R-squared: 0.8766
F-statistic: 44.7 on 13 and 67 DF, p-value: < 2.2e-16
> 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,] 0.06795007 0.135900132 0.932049934
[2,] 0.05626816 0.112536323 0.943731838
[3,] 0.02447001 0.048940016 0.975529992
[4,] 0.52875630 0.942487399 0.471243699
[5,] 0.41345931 0.826918628 0.586540686
[6,] 0.90899950 0.182001002 0.091000501
[7,] 0.96541992 0.069160156 0.034580078
[8,] 0.99320424 0.013591513 0.006795756
[9,] 0.98885529 0.022289411 0.011144705
[10,] 0.98633855 0.027322896 0.013661448
[11,] 0.99067621 0.018647580 0.009323790
[12,] 0.99373713 0.012525747 0.006262874
[13,] 0.99143847 0.017123056 0.008561528
[14,] 0.98665790 0.026684207 0.013342103
[15,] 0.97908235 0.041835300 0.020917650
[16,] 0.98837280 0.023254399 0.011627200
[17,] 0.99143006 0.017139876 0.008569938
[18,] 0.99855100 0.002897993 0.001448996
[19,] 0.99834929 0.003301413 0.001650706
[20,] 0.99741229 0.005175428 0.002587714
[21,] 0.99803056 0.003938887 0.001969443
[22,] 0.99697682 0.006046356 0.003023178
[23,] 0.99681906 0.006361871 0.003180935
[24,] 0.99567469 0.008650629 0.004325315
[25,] 0.99537219 0.009255612 0.004627806
[26,] 0.99325240 0.013495209 0.006747604
[27,] 0.99398111 0.012037788 0.006018894
[28,] 0.99565413 0.008691745 0.004345873
[29,] 0.99600835 0.007983307 0.003991653
[30,] 0.99674148 0.006517049 0.003258525
[31,] 0.99480610 0.010387804 0.005193902
[32,] 0.99253941 0.014921185 0.007460592
[33,] 0.98837746 0.023245077 0.011622539
[34,] 0.98125561 0.037488788 0.018744394
[35,] 0.97697260 0.046054807 0.023027403
[36,] 0.99276218 0.014475633 0.007237817
[37,] 0.99698415 0.006031706 0.003015853
[38,] 0.99594656 0.008106889 0.004053445
[39,] 0.99662041 0.006759185 0.003379593
[40,] 0.99254413 0.014911748 0.007455874
[41,] 0.98679209 0.026415814 0.013207907
[42,] 0.97669655 0.046606894 0.023303447
[43,] 0.96333931 0.073321388 0.036660694
[44,] 0.94531627 0.109367460 0.054683730
[45,] 0.96055028 0.078899450 0.039449725
[46,] 0.96960766 0.060784680 0.030392340
[47,] 0.95302212 0.093955764 0.046977882
[48,] 0.98816133 0.023677332 0.011838666
> postscript(file="/var/www/html/rcomp/tmp/1234b1229175966.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/2vlel1229175966.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/3apnu1229175966.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/4qkkn1229175966.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/54e8d1229175966.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 = 81
Frequency = 1
1 2 3 4 5
7.766609294 4.038037866 3.752323580 4.580895009 10.138037866
6 7 8 9 10
4.723752151 7.998967298 16.970395869 1.098967298 17.473436604
11 12 13 14 15
-3.676563396 -6.476563396 -0.654331612 -0.382903041 -5.668617326
16 17 18 19 20
-2.840045898 1.617096959 -6.597188755 -1.621973609 -5.450545037
21 22 23 24 25
-8.021973609 -11.647504303 -0.697504303 1.202495697 -5.775272519
26 27 28 29 30
-3.703843947 -1.089558233 0.239013196 -9.403843947 -7.318129662
31 32 33 34 35
-7.342914515 2.028514056 -0.642914515 3.931554791 -1.418445209
36 37 38 39 40
-3.618445209 -0.696213425 -3.624784854 -1.610499139 -3.681927711
41 42 43 44 45
-10.824784854 -5.639070568 0.936144578 -1.692426850 0.636144578
46 47 48 49 50
-10.189386116 -1.139386116 -1.139386116 -1.517154332 -1.545725760
51 52 53 54 55
0.868559954 -8.802868617 4.954274240 -0.060011474 -2.384796328
56 57 58 59 60
-2.613367757 3.615203672 3.889672978 6.939672978 8.539672978
61 62 63 64 65
6.261904762 1.133333333 6.047619048 4.276190476 4.133333333
66 67 68 69 70
11.419047619 -1.732243259 -7.560814687 0.167756741 -3.457773953
71 72 73 74 75
-0.007773953 1.492226047 -5.385542169 4.085886403 -2.299827883
76 77 78 79 80
6.228743546 -0.614113597 3.471600688 4.146815835 -1.681755594
81
3.146815835
> postscript(file="/var/www/html/rcomp/tmp/634ku1229175966.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 = 81
Frequency = 1
lag(myerror, k = 1) myerror
0 7.766609294 NA
1 4.038037866 7.766609294
2 3.752323580 4.038037866
3 4.580895009 3.752323580
4 10.138037866 4.580895009
5 4.723752151 10.138037866
6 7.998967298 4.723752151
7 16.970395869 7.998967298
8 1.098967298 16.970395869
9 17.473436604 1.098967298
10 -3.676563396 17.473436604
11 -6.476563396 -3.676563396
12 -0.654331612 -6.476563396
13 -0.382903041 -0.654331612
14 -5.668617326 -0.382903041
15 -2.840045898 -5.668617326
16 1.617096959 -2.840045898
17 -6.597188755 1.617096959
18 -1.621973609 -6.597188755
19 -5.450545037 -1.621973609
20 -8.021973609 -5.450545037
21 -11.647504303 -8.021973609
22 -0.697504303 -11.647504303
23 1.202495697 -0.697504303
24 -5.775272519 1.202495697
25 -3.703843947 -5.775272519
26 -1.089558233 -3.703843947
27 0.239013196 -1.089558233
28 -9.403843947 0.239013196
29 -7.318129662 -9.403843947
30 -7.342914515 -7.318129662
31 2.028514056 -7.342914515
32 -0.642914515 2.028514056
33 3.931554791 -0.642914515
34 -1.418445209 3.931554791
35 -3.618445209 -1.418445209
36 -0.696213425 -3.618445209
37 -3.624784854 -0.696213425
38 -1.610499139 -3.624784854
39 -3.681927711 -1.610499139
40 -10.824784854 -3.681927711
41 -5.639070568 -10.824784854
42 0.936144578 -5.639070568
43 -1.692426850 0.936144578
44 0.636144578 -1.692426850
45 -10.189386116 0.636144578
46 -1.139386116 -10.189386116
47 -1.139386116 -1.139386116
48 -1.517154332 -1.139386116
49 -1.545725760 -1.517154332
50 0.868559954 -1.545725760
51 -8.802868617 0.868559954
52 4.954274240 -8.802868617
53 -0.060011474 4.954274240
54 -2.384796328 -0.060011474
55 -2.613367757 -2.384796328
56 3.615203672 -2.613367757
57 3.889672978 3.615203672
58 6.939672978 3.889672978
59 8.539672978 6.939672978
60 6.261904762 8.539672978
61 1.133333333 6.261904762
62 6.047619048 1.133333333
63 4.276190476 6.047619048
64 4.133333333 4.276190476
65 11.419047619 4.133333333
66 -1.732243259 11.419047619
67 -7.560814687 -1.732243259
68 0.167756741 -7.560814687
69 -3.457773953 0.167756741
70 -0.007773953 -3.457773953
71 1.492226047 -0.007773953
72 -5.385542169 1.492226047
73 4.085886403 -5.385542169
74 -2.299827883 4.085886403
75 6.228743546 -2.299827883
76 -0.614113597 6.228743546
77 3.471600688 -0.614113597
78 4.146815835 3.471600688
79 -1.681755594 4.146815835
80 3.146815835 -1.681755594
81 NA 3.146815835
> dum1 <- dum[2:length(myerror),]
> dum1
lag(myerror, k = 1) myerror
[1,] 4.038037866 7.766609294
[2,] 3.752323580 4.038037866
[3,] 4.580895009 3.752323580
[4,] 10.138037866 4.580895009
[5,] 4.723752151 10.138037866
[6,] 7.998967298 4.723752151
[7,] 16.970395869 7.998967298
[8,] 1.098967298 16.970395869
[9,] 17.473436604 1.098967298
[10,] -3.676563396 17.473436604
[11,] -6.476563396 -3.676563396
[12,] -0.654331612 -6.476563396
[13,] -0.382903041 -0.654331612
[14,] -5.668617326 -0.382903041
[15,] -2.840045898 -5.668617326
[16,] 1.617096959 -2.840045898
[17,] -6.597188755 1.617096959
[18,] -1.621973609 -6.597188755
[19,] -5.450545037 -1.621973609
[20,] -8.021973609 -5.450545037
[21,] -11.647504303 -8.021973609
[22,] -0.697504303 -11.647504303
[23,] 1.202495697 -0.697504303
[24,] -5.775272519 1.202495697
[25,] -3.703843947 -5.775272519
[26,] -1.089558233 -3.703843947
[27,] 0.239013196 -1.089558233
[28,] -9.403843947 0.239013196
[29,] -7.318129662 -9.403843947
[30,] -7.342914515 -7.318129662
[31,] 2.028514056 -7.342914515
[32,] -0.642914515 2.028514056
[33,] 3.931554791 -0.642914515
[34,] -1.418445209 3.931554791
[35,] -3.618445209 -1.418445209
[36,] -0.696213425 -3.618445209
[37,] -3.624784854 -0.696213425
[38,] -1.610499139 -3.624784854
[39,] -3.681927711 -1.610499139
[40,] -10.824784854 -3.681927711
[41,] -5.639070568 -10.824784854
[42,] 0.936144578 -5.639070568
[43,] -1.692426850 0.936144578
[44,] 0.636144578 -1.692426850
[45,] -10.189386116 0.636144578
[46,] -1.139386116 -10.189386116
[47,] -1.139386116 -1.139386116
[48,] -1.517154332 -1.139386116
[49,] -1.545725760 -1.517154332
[50,] 0.868559954 -1.545725760
[51,] -8.802868617 0.868559954
[52,] 4.954274240 -8.802868617
[53,] -0.060011474 4.954274240
[54,] -2.384796328 -0.060011474
[55,] -2.613367757 -2.384796328
[56,] 3.615203672 -2.613367757
[57,] 3.889672978 3.615203672
[58,] 6.939672978 3.889672978
[59,] 8.539672978 6.939672978
[60,] 6.261904762 8.539672978
[61,] 1.133333333 6.261904762
[62,] 6.047619048 1.133333333
[63,] 4.276190476 6.047619048
[64,] 4.133333333 4.276190476
[65,] 11.419047619 4.133333333
[66,] -1.732243259 11.419047619
[67,] -7.560814687 -1.732243259
[68,] 0.167756741 -7.560814687
[69,] -3.457773953 0.167756741
[70,] -0.007773953 -3.457773953
[71,] 1.492226047 -0.007773953
[72,] -5.385542169 1.492226047
[73,] 4.085886403 -5.385542169
[74,] -2.299827883 4.085886403
[75,] 6.228743546 -2.299827883
[76,] -0.614113597 6.228743546
[77,] 3.471600688 -0.614113597
[78,] 4.146815835 3.471600688
[79,] -1.681755594 4.146815835
[80,] 3.146815835 -1.681755594
> z <- as.data.frame(dum1)
> z
lag(myerror, k = 1) myerror
1 4.038037866 7.766609294
2 3.752323580 4.038037866
3 4.580895009 3.752323580
4 10.138037866 4.580895009
5 4.723752151 10.138037866
6 7.998967298 4.723752151
7 16.970395869 7.998967298
8 1.098967298 16.970395869
9 17.473436604 1.098967298
10 -3.676563396 17.473436604
11 -6.476563396 -3.676563396
12 -0.654331612 -6.476563396
13 -0.382903041 -0.654331612
14 -5.668617326 -0.382903041
15 -2.840045898 -5.668617326
16 1.617096959 -2.840045898
17 -6.597188755 1.617096959
18 -1.621973609 -6.597188755
19 -5.450545037 -1.621973609
20 -8.021973609 -5.450545037
21 -11.647504303 -8.021973609
22 -0.697504303 -11.647504303
23 1.202495697 -0.697504303
24 -5.775272519 1.202495697
25 -3.703843947 -5.775272519
26 -1.089558233 -3.703843947
27 0.239013196 -1.089558233
28 -9.403843947 0.239013196
29 -7.318129662 -9.403843947
30 -7.342914515 -7.318129662
31 2.028514056 -7.342914515
32 -0.642914515 2.028514056
33 3.931554791 -0.642914515
34 -1.418445209 3.931554791
35 -3.618445209 -1.418445209
36 -0.696213425 -3.618445209
37 -3.624784854 -0.696213425
38 -1.610499139 -3.624784854
39 -3.681927711 -1.610499139
40 -10.824784854 -3.681927711
41 -5.639070568 -10.824784854
42 0.936144578 -5.639070568
43 -1.692426850 0.936144578
44 0.636144578 -1.692426850
45 -10.189386116 0.636144578
46 -1.139386116 -10.189386116
47 -1.139386116 -1.139386116
48 -1.517154332 -1.139386116
49 -1.545725760 -1.517154332
50 0.868559954 -1.545725760
51 -8.802868617 0.868559954
52 4.954274240 -8.802868617
53 -0.060011474 4.954274240
54 -2.384796328 -0.060011474
55 -2.613367757 -2.384796328
56 3.615203672 -2.613367757
57 3.889672978 3.615203672
58 6.939672978 3.889672978
59 8.539672978 6.939672978
60 6.261904762 8.539672978
61 1.133333333 6.261904762
62 6.047619048 1.133333333
63 4.276190476 6.047619048
64 4.133333333 4.276190476
65 11.419047619 4.133333333
66 -1.732243259 11.419047619
67 -7.560814687 -1.732243259
68 0.167756741 -7.560814687
69 -3.457773953 0.167756741
70 -0.007773953 -3.457773953
71 1.492226047 -0.007773953
72 -5.385542169 1.492226047
73 4.085886403 -5.385542169
74 -2.299827883 4.085886403
75 6.228743546 -2.299827883
76 -0.614113597 6.228743546
77 3.471600688 -0.614113597
78 4.146815835 3.471600688
79 -1.681755594 4.146815835
80 3.146815835 -1.681755594
> 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/7r4o21229175966.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/8jmz51229175966.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/94h991229175966.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/10ygye1229175966.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/11hdsx1229175966.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/126t841229175966.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/13fr301229175966.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/14gaft1229175966.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/15mryd1229175966.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/167sdn1229175967.tab")
+ }
>
> system("convert tmp/1234b1229175966.ps tmp/1234b1229175966.png")
> system("convert tmp/2vlel1229175966.ps tmp/2vlel1229175966.png")
> system("convert tmp/3apnu1229175966.ps tmp/3apnu1229175966.png")
> system("convert tmp/4qkkn1229175966.ps tmp/4qkkn1229175966.png")
> system("convert tmp/54e8d1229175966.ps tmp/54e8d1229175966.png")
> system("convert tmp/634ku1229175966.ps tmp/634ku1229175966.png")
> system("convert tmp/7r4o21229175966.ps tmp/7r4o21229175966.png")
> system("convert tmp/8jmz51229175966.ps tmp/8jmz51229175966.png")
> system("convert tmp/94h991229175966.ps tmp/94h991229175966.png")
> system("convert tmp/10ygye1229175966.ps tmp/10ygye1229175966.png")
>
>
> proc.time()
user system elapsed
9.083 4.670 9.632