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 = '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
Y X M1 M2 M3 M4 M5 M6 M7 M8 M9 M10 M11
1 8.2 20.3 1 0 0 0 0 0 0 0 0 0 0
2 8.0 20.3 0 1 0 0 0 0 0 0 0 0 0
3 7.5 20.3 0 0 1 0 0 0 0 0 0 0 0
4 6.8 15.8 0 0 0 1 0 0 0 0 0 0 0
5 6.5 15.8 0 0 0 0 1 0 0 0 0 0 0
6 6.6 15.8 0 0 0 0 0 1 0 0 0 0 0
7 7.6 23.2 0 0 0 0 0 0 1 0 0 0 0
8 8.0 23.2 0 0 0 0 0 0 0 1 0 0 0
9 8.1 23.2 0 0 0 0 0 0 0 0 1 0 0
10 7.7 20.9 0 0 0 0 0 0 0 0 0 1 0
11 7.5 20.9 0 0 0 0 0 0 0 0 0 0 1
12 7.6 20.9 0 0 0 0 0 0 0 0 0 0 0
13 7.8 19.8 1 0 0 0 0 0 0 0 0 0 0
14 7.8 19.8 0 1 0 0 0 0 0 0 0 0 0
15 7.8 19.8 0 0 1 0 0 0 0 0 0 0 0
16 7.5 20.6 0 0 0 1 0 0 0 0 0 0 0
17 7.5 20.6 0 0 0 0 1 0 0 0 0 0 0
18 7.1 20.6 0 0 0 0 0 1 0 0 0 0 0
19 7.5 21.1 0 0 0 0 0 0 1 0 0 0 0
20 7.5 21.1 0 0 0 0 0 0 0 1 0 0 0
21 7.6 21.1 0 0 0 0 0 0 0 0 1 0 0
22 7.7 22.4 0 0 0 0 0 0 0 0 0 1 0
23 7.7 22.4 0 0 0 0 0 0 0 0 0 0 1
24 7.9 22.4 0 0 0 0 0 0 0 0 0 0 0
25 8.1 20.5 1 0 0 0 0 0 0 0 0 0 0
26 8.2 20.5 0 1 0 0 0 0 0 0 0 0 0
27 8.2 20.5 0 0 1 0 0 0 0 0 0 0 0
28 8.2 18.4 0 0 0 1 0 0 0 0 0 0 0
29 7.9 18.4 0 0 0 0 1 0 0 0 0 0 0
30 7.3 18.4 0 0 0 0 0 1 0 0 0 0 0
31 6.9 17.6 0 0 0 0 0 0 1 0 0 0 0
32 6.6 17.6 0 0 0 0 0 0 0 1 0 0 0
33 6.7 17.6 0 0 0 0 0 0 0 0 1 0 0
34 6.9 18.5 0 0 0 0 0 0 0 0 0 1 0
35 7.0 18.5 0 0 0 0 0 0 0 0 0 0 1
36 7.1 18.5 0 0 0 0 0 0 0 0 0 0 0
37 7.2 17.3 1 0 0 0 0 0 0 0 0 0 0
38 7.1 17.3 0 1 0 0 0 0 0 0 0 0 0
39 6.9 17.3 0 0 1 0 0 0 0 0 0 0 0
40 7.0 16.2 0 0 0 1 0 0 0 0 0 0 0
41 6.8 16.2 0 0 0 0 1 0 0 0 0 0 0
42 6.4 16.2 0 0 0 0 0 1 0 0 0 0 0
43 6.7 18.5 0 0 0 0 0 0 1 0 0 0 0
44 6.6 18.5 0 0 0 0 0 0 0 1 0 0 0
45 6.4 18.5 0 0 0 0 0 0 0 0 1 0 0
46 6.3 16.3 0 0 0 0 0 0 0 0 0 1 0
47 6.2 16.3 0 0 0 0 0 0 0 0 0 0 1
48 6.5 16.3 0 0 0 0 0 0 0 0 0 0 0
49 6.8 16.8 1 0 0 0 0 0 0 0 0 0 0
50 6.8 16.8 0 1 0 0 0 0 0 0 0 0 0
51 6.4 16.8 0 0 1 0 0 0 0 0 0 0 0
52 6.1 14.8 0 0 0 1 0 0 0 0 0 0 0
53 5.8 14.8 0 0 0 0 1 0 0 0 0 0 0
54 6.1 14.8 0 0 0 0 0 1 0 0 0 0 0
55 7.2 21.4 0 0 0 0 0 0 1 0 0 0 0
56 7.3 21.4 0 0 0 0 0 0 0 1 0 0 0
57 6.9 21.4 0 0 0 0 0 0 0 0 1 0 0
58 6.1 16.1 0 0 0 0 0 0 0 0 0 1 0
59 5.8 16.1 0 0 0 0 0 0 0 0 0 0 1
60 6.2 16.1 0 0 0 0 0 0 0 0 0 0 0
61 7.1 19.6 1 0 0 0 0 0 0 0 0 0 0
62 7.7 19.6 0 1 0 0 0 0 0 0 0 0 0
63 7.9 19.6 0 0 1 0 0 0 0 0 0 0 0
64 7.7 18.9 0 0 0 1 0 0 0 0 0 0 0
65 7.4 18.9 0 0 0 0 1 0 0 0 0 0 0
66 7.5 18.9 0 0 0 0 0 1 0 0 0 0 0
67 8.0 21.9 0 0 0 0 0 0 1 0 0 0 0
68 8.1 21.9 0 0 0 0 0 0 0 1 0 0 0
69 8.0 21.9 0 0 0 0 0 0 0 0 1 0 0
> k <- length(x[1,])
> df <- as.data.frame(x)
> (mylm <- lm(df))
Call:
lm(formula = df)
Coefficients:
(Intercept) X M1 M2 M3 M4
1.9294 0.2723 0.4161 0.4828 0.3328 0.5352
M5 M6 M7 M8 M9 M10
0.3019 0.1519 -0.2272 -0.1938 -0.2605 -0.1200
M11
-0.2200
> (mysum <- summary(mylm))
Call:
lm(formula = df)
Residuals:
Min 1Q Median 3Q Max
-0.59665 -0.12947 0.03267 0.17180 0.72463
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 1.9295 0.3783 5.101 4.17e-06 ***
X 0.2723 0.0187 14.565 < 2e-16 ***
M1 0.4162 0.1867 2.229 0.02987 *
M2 0.4828 0.1867 2.586 0.01235 *
M3 0.3328 0.1867 1.782 0.08011 .
M4 0.5352 0.1885 2.839 0.00629 **
M5 0.3019 0.1885 1.602 0.11489
M6 0.1519 0.1885 0.806 0.42383
M7 -0.2272 0.1896 -1.198 0.23597
M8 -0.1938 0.1896 -1.022 0.31108
M9 -0.2605 0.1896 -1.374 0.17498
M10 -0.1200 0.1950 -0.615 0.54076
M11 -0.2200 0.1950 -1.128 0.26401
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 0.3083 on 56 degrees of freedom
Multiple R-squared: 0.822, Adjusted R-squared: 0.7839
F-statistic: 21.56 on 12 and 56 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.252397425 0.504794850 0.7476026
[2,] 0.167602367 0.335204735 0.8323976
[3,] 0.160427153 0.320854307 0.8395728
[4,] 0.107889409 0.215778818 0.8921106
[5,] 0.067332744 0.134665489 0.9326673
[6,] 0.041664671 0.083329342 0.9583353
[7,] 0.031581644 0.063163289 0.9684184
[8,] 0.017025220 0.034050440 0.9829748
[9,] 0.010289976 0.020579953 0.9897100
[10,] 0.004905119 0.009810237 0.9950949
[11,] 0.003892719 0.007785438 0.9961073
[12,] 0.013287698 0.026575395 0.9867123
[13,] 0.332089052 0.664178105 0.6679109
[14,] 0.664902852 0.670194296 0.3350971
[15,] 0.640076689 0.719846621 0.3599233
[16,] 0.641576996 0.716846009 0.3584230
[17,] 0.620531371 0.758937257 0.3794686
[18,] 0.687936206 0.624127588 0.3120638
[19,] 0.620911919 0.758176163 0.3790881
[20,] 0.549463119 0.901073761 0.4505369
[21,] 0.462781866 0.925563732 0.5372181
[22,] 0.512536826 0.974926349 0.4874632
[23,] 0.466625047 0.933250093 0.5333750
[24,] 0.413894086 0.827788173 0.5861059
[25,] 0.369684636 0.739369273 0.6303154
[26,] 0.358403192 0.716806384 0.6415968
[27,] 0.293231784 0.586463568 0.7067682
[28,] 0.243622694 0.487245387 0.7563773
[29,] 0.193600979 0.387201958 0.8063990
[30,] 0.194039123 0.388078247 0.8059609
[31,] 0.138828103 0.277656207 0.8611719
[32,] 0.109496207 0.218992414 0.8905038
[33,] 0.073551081 0.147102162 0.9264489
[34,] 0.144617870 0.289235741 0.8553821
[35,] 0.108643847 0.217287694 0.8913562
[36,] 0.103813389 0.207626778 0.8961866
[37,] 0.069796678 0.139593356 0.9302033
[38,] 0.043878620 0.087757241 0.9561214
> postscript(file="/var/www/html/rcomp/tmp/1mpuw1258665288.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/2mel31258665288.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/3mi9s1258665288.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/4b7ny1258665288.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/5grzo1258665288.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 6
0.326263711 0.059597044 -0.290402956 0.032665235 -0.034001432 0.215998568
7 8 9 10 11 12
-0.420166108 -0.053499442 0.113167225 0.199015929 0.099015929 -0.020984071
13 14 15 16 17 18
0.062424893 -0.004241773 0.145758227 -0.574482115 -0.341148781 -0.591148781
19 20 21 22 23 24
0.051710857 0.018377524 0.185044191 -0.209467618 -0.109467618 -0.129467618
25 26 27 28 29 30
0.171799238 0.205132572 0.355132572 0.724627087 0.657960420 0.207960420
31 32 33 34 35 36
0.404839133 0.071505800 0.238172466 0.052589604 0.252589604 0.132589604
37 38 39 40 41 42
0.143230805 -0.023435862 -0.073435862 0.123736289 0.157069622 -0.092930378
43 44 45 46 47 48
-0.040250995 -0.173584329 -0.306917662 0.051698806 0.051698806 0.131698806
49 50 51 52 53 54
-0.120608013 -0.187274680 -0.437274680 -0.395012401 -0.461679067 -0.011679067
55 56 57 58 59 60
-0.329985852 -0.263319186 -0.596652519 -0.093836721 -0.293836721 -0.113836721
61 62 63 64 65 66
-0.583110634 -0.049777300 0.300222700 0.088465905 0.021799238 0.271799238
67 68 69
0.333852966 0.400519632 0.367186299
> postscript(file="/var/www/html/rcomp/tmp/6tsbb1258665288.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.326263711 NA
1 0.059597044 0.326263711
2 -0.290402956 0.059597044
3 0.032665235 -0.290402956
4 -0.034001432 0.032665235
5 0.215998568 -0.034001432
6 -0.420166108 0.215998568
7 -0.053499442 -0.420166108
8 0.113167225 -0.053499442
9 0.199015929 0.113167225
10 0.099015929 0.199015929
11 -0.020984071 0.099015929
12 0.062424893 -0.020984071
13 -0.004241773 0.062424893
14 0.145758227 -0.004241773
15 -0.574482115 0.145758227
16 -0.341148781 -0.574482115
17 -0.591148781 -0.341148781
18 0.051710857 -0.591148781
19 0.018377524 0.051710857
20 0.185044191 0.018377524
21 -0.209467618 0.185044191
22 -0.109467618 -0.209467618
23 -0.129467618 -0.109467618
24 0.171799238 -0.129467618
25 0.205132572 0.171799238
26 0.355132572 0.205132572
27 0.724627087 0.355132572
28 0.657960420 0.724627087
29 0.207960420 0.657960420
30 0.404839133 0.207960420
31 0.071505800 0.404839133
32 0.238172466 0.071505800
33 0.052589604 0.238172466
34 0.252589604 0.052589604
35 0.132589604 0.252589604
36 0.143230805 0.132589604
37 -0.023435862 0.143230805
38 -0.073435862 -0.023435862
39 0.123736289 -0.073435862
40 0.157069622 0.123736289
41 -0.092930378 0.157069622
42 -0.040250995 -0.092930378
43 -0.173584329 -0.040250995
44 -0.306917662 -0.173584329
45 0.051698806 -0.306917662
46 0.051698806 0.051698806
47 0.131698806 0.051698806
48 -0.120608013 0.131698806
49 -0.187274680 -0.120608013
50 -0.437274680 -0.187274680
51 -0.395012401 -0.437274680
52 -0.461679067 -0.395012401
53 -0.011679067 -0.461679067
54 -0.329985852 -0.011679067
55 -0.263319186 -0.329985852
56 -0.596652519 -0.263319186
57 -0.093836721 -0.596652519
58 -0.293836721 -0.093836721
59 -0.113836721 -0.293836721
60 -0.583110634 -0.113836721
61 -0.049777300 -0.583110634
62 0.300222700 -0.049777300
63 0.088465905 0.300222700
64 0.021799238 0.088465905
65 0.271799238 0.021799238
66 0.333852966 0.271799238
67 0.400519632 0.333852966
68 0.367186299 0.400519632
69 NA 0.367186299
> dum1 <- dum[2:length(myerror),]
> dum1
lag(myerror, k = 1) myerror
[1,] 0.059597044 0.326263711
[2,] -0.290402956 0.059597044
[3,] 0.032665235 -0.290402956
[4,] -0.034001432 0.032665235
[5,] 0.215998568 -0.034001432
[6,] -0.420166108 0.215998568
[7,] -0.053499442 -0.420166108
[8,] 0.113167225 -0.053499442
[9,] 0.199015929 0.113167225
[10,] 0.099015929 0.199015929
[11,] -0.020984071 0.099015929
[12,] 0.062424893 -0.020984071
[13,] -0.004241773 0.062424893
[14,] 0.145758227 -0.004241773
[15,] -0.574482115 0.145758227
[16,] -0.341148781 -0.574482115
[17,] -0.591148781 -0.341148781
[18,] 0.051710857 -0.591148781
[19,] 0.018377524 0.051710857
[20,] 0.185044191 0.018377524
[21,] -0.209467618 0.185044191
[22,] -0.109467618 -0.209467618
[23,] -0.129467618 -0.109467618
[24,] 0.171799238 -0.129467618
[25,] 0.205132572 0.171799238
[26,] 0.355132572 0.205132572
[27,] 0.724627087 0.355132572
[28,] 0.657960420 0.724627087
[29,] 0.207960420 0.657960420
[30,] 0.404839133 0.207960420
[31,] 0.071505800 0.404839133
[32,] 0.238172466 0.071505800
[33,] 0.052589604 0.238172466
[34,] 0.252589604 0.052589604
[35,] 0.132589604 0.252589604
[36,] 0.143230805 0.132589604
[37,] -0.023435862 0.143230805
[38,] -0.073435862 -0.023435862
[39,] 0.123736289 -0.073435862
[40,] 0.157069622 0.123736289
[41,] -0.092930378 0.157069622
[42,] -0.040250995 -0.092930378
[43,] -0.173584329 -0.040250995
[44,] -0.306917662 -0.173584329
[45,] 0.051698806 -0.306917662
[46,] 0.051698806 0.051698806
[47,] 0.131698806 0.051698806
[48,] -0.120608013 0.131698806
[49,] -0.187274680 -0.120608013
[50,] -0.437274680 -0.187274680
[51,] -0.395012401 -0.437274680
[52,] -0.461679067 -0.395012401
[53,] -0.011679067 -0.461679067
[54,] -0.329985852 -0.011679067
[55,] -0.263319186 -0.329985852
[56,] -0.596652519 -0.263319186
[57,] -0.093836721 -0.596652519
[58,] -0.293836721 -0.093836721
[59,] -0.113836721 -0.293836721
[60,] -0.583110634 -0.113836721
[61,] -0.049777300 -0.583110634
[62,] 0.300222700 -0.049777300
[63,] 0.088465905 0.300222700
[64,] 0.021799238 0.088465905
[65,] 0.271799238 0.021799238
[66,] 0.333852966 0.271799238
[67,] 0.400519632 0.333852966
[68,] 0.367186299 0.400519632
> z <- as.data.frame(dum1)
> z
lag(myerror, k = 1) myerror
1 0.059597044 0.326263711
2 -0.290402956 0.059597044
3 0.032665235 -0.290402956
4 -0.034001432 0.032665235
5 0.215998568 -0.034001432
6 -0.420166108 0.215998568
7 -0.053499442 -0.420166108
8 0.113167225 -0.053499442
9 0.199015929 0.113167225
10 0.099015929 0.199015929
11 -0.020984071 0.099015929
12 0.062424893 -0.020984071
13 -0.004241773 0.062424893
14 0.145758227 -0.004241773
15 -0.574482115 0.145758227
16 -0.341148781 -0.574482115
17 -0.591148781 -0.341148781
18 0.051710857 -0.591148781
19 0.018377524 0.051710857
20 0.185044191 0.018377524
21 -0.209467618 0.185044191
22 -0.109467618 -0.209467618
23 -0.129467618 -0.109467618
24 0.171799238 -0.129467618
25 0.205132572 0.171799238
26 0.355132572 0.205132572
27 0.724627087 0.355132572
28 0.657960420 0.724627087
29 0.207960420 0.657960420
30 0.404839133 0.207960420
31 0.071505800 0.404839133
32 0.238172466 0.071505800
33 0.052589604 0.238172466
34 0.252589604 0.052589604
35 0.132589604 0.252589604
36 0.143230805 0.132589604
37 -0.023435862 0.143230805
38 -0.073435862 -0.023435862
39 0.123736289 -0.073435862
40 0.157069622 0.123736289
41 -0.092930378 0.157069622
42 -0.040250995 -0.092930378
43 -0.173584329 -0.040250995
44 -0.306917662 -0.173584329
45 0.051698806 -0.306917662
46 0.051698806 0.051698806
47 0.131698806 0.051698806
48 -0.120608013 0.131698806
49 -0.187274680 -0.120608013
50 -0.437274680 -0.187274680
51 -0.395012401 -0.437274680
52 -0.461679067 -0.395012401
53 -0.011679067 -0.461679067
54 -0.329985852 -0.011679067
55 -0.263319186 -0.329985852
56 -0.596652519 -0.263319186
57 -0.093836721 -0.596652519
58 -0.293836721 -0.093836721
59 -0.113836721 -0.293836721
60 -0.583110634 -0.113836721
61 -0.049777300 -0.583110634
62 0.300222700 -0.049777300
63 0.088465905 0.300222700
64 0.021799238 0.088465905
65 0.271799238 0.021799238
66 0.333852966 0.271799238
67 0.400519632 0.333852966
68 0.367186299 0.400519632
> 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/76exp1258665288.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/8cxi01258665288.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/9bn2f1258665288.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/10a9yp1258665288.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/11xdpe1258665288.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/12q23c1258665288.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/13jvf31258665288.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/14w4bg1258665288.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/157rm41258665288.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/16kj1d1258665288.tab")
+ }
>
> system("convert tmp/1mpuw1258665288.ps tmp/1mpuw1258665288.png")
> system("convert tmp/2mel31258665288.ps tmp/2mel31258665288.png")
> system("convert tmp/3mi9s1258665288.ps tmp/3mi9s1258665288.png")
> system("convert tmp/4b7ny1258665288.ps tmp/4b7ny1258665288.png")
> system("convert tmp/5grzo1258665288.ps tmp/5grzo1258665288.png")
> system("convert tmp/6tsbb1258665288.ps tmp/6tsbb1258665288.png")
> system("convert tmp/76exp1258665288.ps tmp/76exp1258665288.png")
> system("convert tmp/8cxi01258665288.ps tmp/8cxi01258665288.png")
> system("convert tmp/9bn2f1258665288.ps tmp/9bn2f1258665288.png")
> system("convert tmp/10a9yp1258665288.ps tmp/10a9yp1258665288.png")
>
>
> proc.time()
user system elapsed
2.539 1.567 3.046