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(8.2,20.3,8,20.3,7.5,20.3,6.8,15.8,6.5,15.8,6.6,15.8,7.6,23.2,8,23.2,8.1,23.2,7.7,20.9,7.5,20.9,7.6,20.9,7.8,19.8,7.8,19.8,7.8,19.8,7.5,20.6,7.5,20.6,7.1,20.6,7.5,21.1,7.5,21.1,7.6,21.1,7.7,22.4,7.7,22.4,7.9,22.4,8.1,20.5,8.2,20.5,8.2,20.5,8.2,18.4,7.9,18.4,7.3,18.4,6.9,17.6,6.6,17.6,6.7,17.6,6.9,18.5,7,18.5,7.1,18.5,7.2,17.3,7.1,17.3,6.9,17.3,7,16.2,6.8,16.2,6.4,16.2,6.7,18.5,6.6,18.5,6.4,18.5,6.3,16.3,6.2,16.3,6.5,16.3,6.8,16.8,6.8,16.8,6.4,16.8,6.1,14.8,5.8,14.8,6.1,14.8,7.2,21.4,7.3,21.4,6.9,21.4,6.1,16.1,5.8,16.1,6.2,16.1,7.1,19.6,7.7,19.6,7.9,19.6,7.7,18.9,7.4,18.9,7.5,18.9,8,21.9,8.1,21.9,8,21.9),dim=c(2,69),dimnames=list(c('Y','X'),1:69))
> y <- array(NA,dim=c(2,69),dimnames=list(c('Y','X'),1:69))
> 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 8.2 20.3 1 0 0 0 0 0 0 0 0 0 0 1
2 8.0 20.3 0 1 0 0 0 0 0 0 0 0 0 2
3 7.5 20.3 0 0 1 0 0 0 0 0 0 0 0 3
4 6.8 15.8 0 0 0 1 0 0 0 0 0 0 0 4
5 6.5 15.8 0 0 0 0 1 0 0 0 0 0 0 5
6 6.6 15.8 0 0 0 0 0 1 0 0 0 0 0 6
7 7.6 23.2 0 0 0 0 0 0 1 0 0 0 0 7
8 8.0 23.2 0 0 0 0 0 0 0 1 0 0 0 8
9 8.1 23.2 0 0 0 0 0 0 0 0 1 0 0 9
10 7.7 20.9 0 0 0 0 0 0 0 0 0 1 0 10
11 7.5 20.9 0 0 0 0 0 0 0 0 0 0 1 11
12 7.6 20.9 0 0 0 0 0 0 0 0 0 0 0 12
13 7.8 19.8 1 0 0 0 0 0 0 0 0 0 0 13
14 7.8 19.8 0 1 0 0 0 0 0 0 0 0 0 14
15 7.8 19.8 0 0 1 0 0 0 0 0 0 0 0 15
16 7.5 20.6 0 0 0 1 0 0 0 0 0 0 0 16
17 7.5 20.6 0 0 0 0 1 0 0 0 0 0 0 17
18 7.1 20.6 0 0 0 0 0 1 0 0 0 0 0 18
19 7.5 21.1 0 0 0 0 0 0 1 0 0 0 0 19
20 7.5 21.1 0 0 0 0 0 0 0 1 0 0 0 20
21 7.6 21.1 0 0 0 0 0 0 0 0 1 0 0 21
22 7.7 22.4 0 0 0 0 0 0 0 0 0 1 0 22
23 7.7 22.4 0 0 0 0 0 0 0 0 0 0 1 23
24 7.9 22.4 0 0 0 0 0 0 0 0 0 0 0 24
25 8.1 20.5 1 0 0 0 0 0 0 0 0 0 0 25
26 8.2 20.5 0 1 0 0 0 0 0 0 0 0 0 26
27 8.2 20.5 0 0 1 0 0 0 0 0 0 0 0 27
28 8.2 18.4 0 0 0 1 0 0 0 0 0 0 0 28
29 7.9 18.4 0 0 0 0 1 0 0 0 0 0 0 29
30 7.3 18.4 0 0 0 0 0 1 0 0 0 0 0 30
31 6.9 17.6 0 0 0 0 0 0 1 0 0 0 0 31
32 6.6 17.6 0 0 0 0 0 0 0 1 0 0 0 32
33 6.7 17.6 0 0 0 0 0 0 0 0 1 0 0 33
34 6.9 18.5 0 0 0 0 0 0 0 0 0 1 0 34
35 7.0 18.5 0 0 0 0 0 0 0 0 0 0 1 35
36 7.1 18.5 0 0 0 0 0 0 0 0 0 0 0 36
37 7.2 17.3 1 0 0 0 0 0 0 0 0 0 0 37
38 7.1 17.3 0 1 0 0 0 0 0 0 0 0 0 38
39 6.9 17.3 0 0 1 0 0 0 0 0 0 0 0 39
40 7.0 16.2 0 0 0 1 0 0 0 0 0 0 0 40
41 6.8 16.2 0 0 0 0 1 0 0 0 0 0 0 41
42 6.4 16.2 0 0 0 0 0 1 0 0 0 0 0 42
43 6.7 18.5 0 0 0 0 0 0 1 0 0 0 0 43
44 6.6 18.5 0 0 0 0 0 0 0 1 0 0 0 44
45 6.4 18.5 0 0 0 0 0 0 0 0 1 0 0 45
46 6.3 16.3 0 0 0 0 0 0 0 0 0 1 0 46
47 6.2 16.3 0 0 0 0 0 0 0 0 0 0 1 47
48 6.5 16.3 0 0 0 0 0 0 0 0 0 0 0 48
49 6.8 16.8 1 0 0 0 0 0 0 0 0 0 0 49
50 6.8 16.8 0 1 0 0 0 0 0 0 0 0 0 50
51 6.4 16.8 0 0 1 0 0 0 0 0 0 0 0 51
52 6.1 14.8 0 0 0 1 0 0 0 0 0 0 0 52
53 5.8 14.8 0 0 0 0 1 0 0 0 0 0 0 53
54 6.1 14.8 0 0 0 0 0 1 0 0 0 0 0 54
55 7.2 21.4 0 0 0 0 0 0 1 0 0 0 0 55
56 7.3 21.4 0 0 0 0 0 0 0 1 0 0 0 56
57 6.9 21.4 0 0 0 0 0 0 0 0 1 0 0 57
58 6.1 16.1 0 0 0 0 0 0 0 0 0 1 0 58
59 5.8 16.1 0 0 0 0 0 0 0 0 0 0 1 59
60 6.2 16.1 0 0 0 0 0 0 0 0 0 0 0 60
61 7.1 19.6 1 0 0 0 0 0 0 0 0 0 0 61
62 7.7 19.6 0 1 0 0 0 0 0 0 0 0 0 62
63 7.9 19.6 0 0 1 0 0 0 0 0 0 0 0 63
64 7.7 18.9 0 0 0 1 0 0 0 0 0 0 0 64
65 7.4 18.9 0 0 0 0 1 0 0 0 0 0 0 65
66 7.5 18.9 0 0 0 0 0 1 0 0 0 0 0 66
67 8.0 21.9 0 0 0 0 0 0 1 0 0 0 0 67
68 8.1 21.9 0 0 0 0 0 0 0 1 0 0 0 68
69 8.0 21.9 0 0 0 0 0 0 0 0 1 0 0 69
> k <- length(x[1,])
> df <- as.data.frame(x)
> (mylm <- lm(df))
Call:
lm(formula = df)
Coefficients:
(Intercept) X M1 M2 M3 M4
2.047491 0.268190 0.411431 0.479214 0.330330 0.527218
M5 M6 M7 M8 M9 M10
0.295001 0.146118 -0.218702 -0.184252 -0.249802 -0.122233
M11 t
-0.221117 -0.001117
> (mysum <- summary(mylm))
Call:
lm(formula = df)
Residuals:
Min 1Q Median 3Q Max
-0.59823 -0.12816 0.03127 0.16694 0.72185
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 2.047491 0.437171 4.684 1.89e-05 ***
X 0.268190 0.020264 13.235 < 2e-16 ***
M1 0.411431 0.188097 2.187 0.0330 *
M2 0.479214 0.188015 2.549 0.0136 *
M3 0.330330 0.187956 1.757 0.0844 .
M4 0.527218 0.190227 2.772 0.0076 **
M5 0.295001 0.190082 1.552 0.1264
M6 0.146118 0.189960 0.769 0.4451
M7 -0.218702 0.191431 -1.142 0.2582
M8 -0.184252 0.191606 -0.962 0.3404
M9 -0.249802 0.191801 -1.302 0.1982
M10 -0.122233 0.196255 -0.623 0.5360
M11 -0.221117 0.196224 -1.127 0.2647
t -0.001117 0.002033 -0.549 0.5852
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 0.3102 on 55 degrees of freedom
Multiple R-squared: 0.823, Adjusted R-squared: 0.7812
F-statistic: 19.67 on 13 and 55 DF, p-value: 3.908e-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.30565191 0.61130382 0.69434809
[2,] 0.28540861 0.57081723 0.71459139
[3,] 0.19666347 0.39332694 0.80333653
[4,] 0.13487222 0.26974444 0.86512778
[5,] 0.08237132 0.16474263 0.91762868
[6,] 0.06327432 0.12654864 0.93672568
[7,] 0.04048949 0.08097899 0.95951051
[8,] 0.03688079 0.07376157 0.96311921
[9,] 0.02009080 0.04018160 0.97990920
[10,] 0.01780545 0.03561091 0.98219455
[11,] 0.03688346 0.07376693 0.96311654
[12,] 0.32347519 0.64695037 0.67652481
[13,] 0.48051549 0.96103098 0.51948451
[14,] 0.40065061 0.80130121 0.59934939
[15,] 0.41032280 0.82064560 0.58967720
[16,] 0.49001760 0.98003519 0.50998240
[17,] 0.58706548 0.82586904 0.41293452
[18,] 0.54920425 0.90159149 0.45079575
[19,] 0.46552007 0.93104014 0.53447993
[20,] 0.39544417 0.79088835 0.60455583
[21,] 0.44899312 0.89798623 0.55100688
[22,] 0.42119356 0.84238713 0.57880644
[23,] 0.38111217 0.76222435 0.61888783
[24,] 0.33988516 0.67977031 0.66011484
[25,] 0.36319056 0.72638112 0.63680944
[26,] 0.28806072 0.57612144 0.71193928
[27,] 0.24877910 0.49755819 0.75122090
[28,] 0.20110782 0.40221564 0.79889218
[29,] 0.20589221 0.41178442 0.79410779
[30,] 0.16782761 0.33565523 0.83217239
[31,] 0.22035505 0.44071009 0.77964495
[32,] 0.41394157 0.82788315 0.58605843
[33,] 0.92996031 0.14007938 0.07003969
[34,] 0.97721474 0.04557051 0.02278526
[35,] 0.96817173 0.06365654 0.03182827
[36,] 0.92134207 0.15731585 0.07865793
> postscript(file="/var/www/html/rcomp/tmp/1fjeg1258666947.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/2bpjw1258666947.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/3bzhj1258666947.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/4rbdg1258666947.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/5hfcm1258666947.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 = 69
Frequency = 1
1 2 3 4 5
0.2979322160 0.0312655494 -0.3187344506 -0.0076492978 -0.0743159645
6 7 8 9 10
0.1756840355 -0.4429881608 -0.0763214941 0.0903451726 0.1807307399
11 12 13 14 15
0.0807307399 -0.0392692601 0.0454259963 -0.0212406704 0.1287593296
16 17 18 19 20
-0.5815640153 -0.3482306819 -0.5982306819 0.0336100716 0.0002767383
21 22 23 24 25
0.1669434050 -0.2081560450 -0.1081560450 -0.1281560450 0.1710914374
26 27 28 29 30
0.2044247708 0.3544247708 0.7218532454 0.6551865787 0.2051865787
31 32 33 34 35
0.3856746997 0.0523413663 0.2190080330 0.0511846961 0.2511846961
36 37 38 39 40
0.1311846961 0.1426989807 -0.0239676860 -0.0739676860 0.1252705060
41 42 43 44 45
0.1586038394 -0.0913961606 -0.0422979157 -0.1756312490 -0.3089645824
46 47 48 49 50
0.0546019567 0.0546019567 0.1346019567 -0.1098072391 -0.1764739057
51 52 53 54 55
-0.4264739057 -0.3858644594 -0.4525311260 -0.0025311260 -0.3066510963
56 57 58 59 60
-0.2399844296 -0.5733177629 -0.0783613478 -0.2783613478 -0.0983613478
61 62 63 64 65
-0.5473413914 -0.0140080580 0.3359919420 0.1279540210 0.0612873543
66 67 68 69
0.3112873543 0.3726524014 0.4393190681 0.4059857347
> postscript(file="/var/www/html/rcomp/tmp/613ue1258666947.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 = 69
Frequency = 1
lag(myerror, k = 1) myerror
0 0.2979322160 NA
1 0.0312655494 0.2979322160
2 -0.3187344506 0.0312655494
3 -0.0076492978 -0.3187344506
4 -0.0743159645 -0.0076492978
5 0.1756840355 -0.0743159645
6 -0.4429881608 0.1756840355
7 -0.0763214941 -0.4429881608
8 0.0903451726 -0.0763214941
9 0.1807307399 0.0903451726
10 0.0807307399 0.1807307399
11 -0.0392692601 0.0807307399
12 0.0454259963 -0.0392692601
13 -0.0212406704 0.0454259963
14 0.1287593296 -0.0212406704
15 -0.5815640153 0.1287593296
16 -0.3482306819 -0.5815640153
17 -0.5982306819 -0.3482306819
18 0.0336100716 -0.5982306819
19 0.0002767383 0.0336100716
20 0.1669434050 0.0002767383
21 -0.2081560450 0.1669434050
22 -0.1081560450 -0.2081560450
23 -0.1281560450 -0.1081560450
24 0.1710914374 -0.1281560450
25 0.2044247708 0.1710914374
26 0.3544247708 0.2044247708
27 0.7218532454 0.3544247708
28 0.6551865787 0.7218532454
29 0.2051865787 0.6551865787
30 0.3856746997 0.2051865787
31 0.0523413663 0.3856746997
32 0.2190080330 0.0523413663
33 0.0511846961 0.2190080330
34 0.2511846961 0.0511846961
35 0.1311846961 0.2511846961
36 0.1426989807 0.1311846961
37 -0.0239676860 0.1426989807
38 -0.0739676860 -0.0239676860
39 0.1252705060 -0.0739676860
40 0.1586038394 0.1252705060
41 -0.0913961606 0.1586038394
42 -0.0422979157 -0.0913961606
43 -0.1756312490 -0.0422979157
44 -0.3089645824 -0.1756312490
45 0.0546019567 -0.3089645824
46 0.0546019567 0.0546019567
47 0.1346019567 0.0546019567
48 -0.1098072391 0.1346019567
49 -0.1764739057 -0.1098072391
50 -0.4264739057 -0.1764739057
51 -0.3858644594 -0.4264739057
52 -0.4525311260 -0.3858644594
53 -0.0025311260 -0.4525311260
54 -0.3066510963 -0.0025311260
55 -0.2399844296 -0.3066510963
56 -0.5733177629 -0.2399844296
57 -0.0783613478 -0.5733177629
58 -0.2783613478 -0.0783613478
59 -0.0983613478 -0.2783613478
60 -0.5473413914 -0.0983613478
61 -0.0140080580 -0.5473413914
62 0.3359919420 -0.0140080580
63 0.1279540210 0.3359919420
64 0.0612873543 0.1279540210
65 0.3112873543 0.0612873543
66 0.3726524014 0.3112873543
67 0.4393190681 0.3726524014
68 0.4059857347 0.4393190681
69 NA 0.4059857347
> dum1 <- dum[2:length(myerror),]
> dum1
lag(myerror, k = 1) myerror
[1,] 0.0312655494 0.2979322160
[2,] -0.3187344506 0.0312655494
[3,] -0.0076492978 -0.3187344506
[4,] -0.0743159645 -0.0076492978
[5,] 0.1756840355 -0.0743159645
[6,] -0.4429881608 0.1756840355
[7,] -0.0763214941 -0.4429881608
[8,] 0.0903451726 -0.0763214941
[9,] 0.1807307399 0.0903451726
[10,] 0.0807307399 0.1807307399
[11,] -0.0392692601 0.0807307399
[12,] 0.0454259963 -0.0392692601
[13,] -0.0212406704 0.0454259963
[14,] 0.1287593296 -0.0212406704
[15,] -0.5815640153 0.1287593296
[16,] -0.3482306819 -0.5815640153
[17,] -0.5982306819 -0.3482306819
[18,] 0.0336100716 -0.5982306819
[19,] 0.0002767383 0.0336100716
[20,] 0.1669434050 0.0002767383
[21,] -0.2081560450 0.1669434050
[22,] -0.1081560450 -0.2081560450
[23,] -0.1281560450 -0.1081560450
[24,] 0.1710914374 -0.1281560450
[25,] 0.2044247708 0.1710914374
[26,] 0.3544247708 0.2044247708
[27,] 0.7218532454 0.3544247708
[28,] 0.6551865787 0.7218532454
[29,] 0.2051865787 0.6551865787
[30,] 0.3856746997 0.2051865787
[31,] 0.0523413663 0.3856746997
[32,] 0.2190080330 0.0523413663
[33,] 0.0511846961 0.2190080330
[34,] 0.2511846961 0.0511846961
[35,] 0.1311846961 0.2511846961
[36,] 0.1426989807 0.1311846961
[37,] -0.0239676860 0.1426989807
[38,] -0.0739676860 -0.0239676860
[39,] 0.1252705060 -0.0739676860
[40,] 0.1586038394 0.1252705060
[41,] -0.0913961606 0.1586038394
[42,] -0.0422979157 -0.0913961606
[43,] -0.1756312490 -0.0422979157
[44,] -0.3089645824 -0.1756312490
[45,] 0.0546019567 -0.3089645824
[46,] 0.0546019567 0.0546019567
[47,] 0.1346019567 0.0546019567
[48,] -0.1098072391 0.1346019567
[49,] -0.1764739057 -0.1098072391
[50,] -0.4264739057 -0.1764739057
[51,] -0.3858644594 -0.4264739057
[52,] -0.4525311260 -0.3858644594
[53,] -0.0025311260 -0.4525311260
[54,] -0.3066510963 -0.0025311260
[55,] -0.2399844296 -0.3066510963
[56,] -0.5733177629 -0.2399844296
[57,] -0.0783613478 -0.5733177629
[58,] -0.2783613478 -0.0783613478
[59,] -0.0983613478 -0.2783613478
[60,] -0.5473413914 -0.0983613478
[61,] -0.0140080580 -0.5473413914
[62,] 0.3359919420 -0.0140080580
[63,] 0.1279540210 0.3359919420
[64,] 0.0612873543 0.1279540210
[65,] 0.3112873543 0.0612873543
[66,] 0.3726524014 0.3112873543
[67,] 0.4393190681 0.3726524014
[68,] 0.4059857347 0.4393190681
> z <- as.data.frame(dum1)
> z
lag(myerror, k = 1) myerror
1 0.0312655494 0.2979322160
2 -0.3187344506 0.0312655494
3 -0.0076492978 -0.3187344506
4 -0.0743159645 -0.0076492978
5 0.1756840355 -0.0743159645
6 -0.4429881608 0.1756840355
7 -0.0763214941 -0.4429881608
8 0.0903451726 -0.0763214941
9 0.1807307399 0.0903451726
10 0.0807307399 0.1807307399
11 -0.0392692601 0.0807307399
12 0.0454259963 -0.0392692601
13 -0.0212406704 0.0454259963
14 0.1287593296 -0.0212406704
15 -0.5815640153 0.1287593296
16 -0.3482306819 -0.5815640153
17 -0.5982306819 -0.3482306819
18 0.0336100716 -0.5982306819
19 0.0002767383 0.0336100716
20 0.1669434050 0.0002767383
21 -0.2081560450 0.1669434050
22 -0.1081560450 -0.2081560450
23 -0.1281560450 -0.1081560450
24 0.1710914374 -0.1281560450
25 0.2044247708 0.1710914374
26 0.3544247708 0.2044247708
27 0.7218532454 0.3544247708
28 0.6551865787 0.7218532454
29 0.2051865787 0.6551865787
30 0.3856746997 0.2051865787
31 0.0523413663 0.3856746997
32 0.2190080330 0.0523413663
33 0.0511846961 0.2190080330
34 0.2511846961 0.0511846961
35 0.1311846961 0.2511846961
36 0.1426989807 0.1311846961
37 -0.0239676860 0.1426989807
38 -0.0739676860 -0.0239676860
39 0.1252705060 -0.0739676860
40 0.1586038394 0.1252705060
41 -0.0913961606 0.1586038394
42 -0.0422979157 -0.0913961606
43 -0.1756312490 -0.0422979157
44 -0.3089645824 -0.1756312490
45 0.0546019567 -0.3089645824
46 0.0546019567 0.0546019567
47 0.1346019567 0.0546019567
48 -0.1098072391 0.1346019567
49 -0.1764739057 -0.1098072391
50 -0.4264739057 -0.1764739057
51 -0.3858644594 -0.4264739057
52 -0.4525311260 -0.3858644594
53 -0.0025311260 -0.4525311260
54 -0.3066510963 -0.0025311260
55 -0.2399844296 -0.3066510963
56 -0.5733177629 -0.2399844296
57 -0.0783613478 -0.5733177629
58 -0.2783613478 -0.0783613478
59 -0.0983613478 -0.2783613478
60 -0.5473413914 -0.0983613478
61 -0.0140080580 -0.5473413914
62 0.3359919420 -0.0140080580
63 0.1279540210 0.3359919420
64 0.0612873543 0.1279540210
65 0.3112873543 0.0612873543
66 0.3726524014 0.3112873543
67 0.4393190681 0.3726524014
68 0.4059857347 0.4393190681
> 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/7iahf1258666947.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/8pmbp1258666947.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/94nf41258666947.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/10xm701258666947.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/11bhec1258666947.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/125v4y1258666947.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/1394851258666947.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/141oli1258666947.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/15v14s1258666947.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/16atwo1258666947.tab")
+ }
>
> system("convert tmp/1fjeg1258666947.ps tmp/1fjeg1258666947.png")
> system("convert tmp/2bpjw1258666947.ps tmp/2bpjw1258666947.png")
> system("convert tmp/3bzhj1258666947.ps tmp/3bzhj1258666947.png")
> system("convert tmp/4rbdg1258666947.ps tmp/4rbdg1258666947.png")
> system("convert tmp/5hfcm1258666947.ps tmp/5hfcm1258666947.png")
> system("convert tmp/613ue1258666947.ps tmp/613ue1258666947.png")
> system("convert tmp/7iahf1258666947.ps tmp/7iahf1258666947.png")
> system("convert tmp/8pmbp1258666947.ps tmp/8pmbp1258666947.png")
> system("convert tmp/94nf41258666947.ps tmp/94nf41258666947.png")
> system("convert tmp/10xm701258666947.ps tmp/10xm701258666947.png")
>
>
> proc.time()
user system elapsed
2.596 1.623 5.915