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(54.8
+ ,0
+ ,56
+ ,56.6
+ ,52.7
+ ,0
+ ,54.8
+ ,56
+ ,50.9
+ ,0
+ ,52.7
+ ,54.8
+ ,50.6
+ ,0
+ ,50.9
+ ,52.7
+ ,52.1
+ ,0
+ ,50.6
+ ,50.9
+ ,53.3
+ ,0
+ ,52.1
+ ,50.6
+ ,53.9
+ ,0
+ ,53.3
+ ,52.1
+ ,54.3
+ ,0
+ ,53.9
+ ,53.3
+ ,54.2
+ ,0
+ ,54.3
+ ,53.9
+ ,54.2
+ ,0
+ ,54.2
+ ,54.3
+ ,53.5
+ ,0
+ ,54.2
+ ,54.2
+ ,51.4
+ ,0
+ ,53.5
+ ,54.2
+ ,50.5
+ ,0
+ ,51.4
+ ,53.5
+ ,50.3
+ ,0
+ ,50.5
+ ,51.4
+ ,49.8
+ ,0
+ ,50.3
+ ,50.5
+ ,50.7
+ ,0
+ ,49.8
+ ,50.3
+ ,52.8
+ ,0
+ ,50.7
+ ,49.8
+ ,55.3
+ ,0
+ ,52.8
+ ,50.7
+ ,57.3
+ ,0
+ ,55.3
+ ,52.8
+ ,57.5
+ ,0
+ ,57.3
+ ,55.3
+ ,56.8
+ ,0
+ ,57.5
+ ,57.3
+ ,56.4
+ ,0
+ ,56.8
+ ,57.5
+ ,56.3
+ ,0
+ ,56.4
+ ,56.8
+ ,56.4
+ ,0
+ ,56.3
+ ,56.4
+ ,57
+ ,0
+ ,56.4
+ ,56.3
+ ,57.9
+ ,0
+ ,57
+ ,56.4
+ ,58.9
+ ,0
+ ,57.9
+ ,57
+ ,58.8
+ ,0
+ ,58.9
+ ,57.9
+ ,56.5
+ ,1
+ ,58.8
+ ,58.9
+ ,51.9
+ ,1
+ ,56.5
+ ,58.8
+ ,47.4
+ ,1
+ ,51.9
+ ,56.5
+ ,44.9
+ ,1
+ ,47.4
+ ,51.9
+ ,43.9
+ ,1
+ ,44.9
+ ,47.4
+ ,43.4
+ ,1
+ ,43.9
+ ,44.9
+ ,42.9
+ ,1
+ ,43.4
+ ,43.9
+ ,42.6
+ ,1
+ ,42.9
+ ,43.4
+ ,42.2
+ ,1
+ ,42.6
+ ,42.9
+ ,41.2
+ ,1
+ ,42.2
+ ,42.6
+ ,40.2
+ ,1
+ ,41.2
+ ,42.2
+ ,39.3
+ ,1
+ ,40.2
+ ,41.2
+ ,38.5
+ ,1
+ ,39.3
+ ,40.2
+ ,38.3
+ ,1
+ ,38.5
+ ,39.3
+ ,37.9
+ ,1
+ ,38.3
+ ,38.5
+ ,37.6
+ ,1
+ ,37.9
+ ,38.3
+ ,37.3
+ ,1
+ ,37.6
+ ,37.9
+ ,36
+ ,1
+ ,37.3
+ ,37.6
+ ,34.5
+ ,1
+ ,36
+ ,37.3
+ ,33.5
+ ,1
+ ,34.5
+ ,36
+ ,32.9
+ ,1
+ ,33.5
+ ,34.5
+ ,32.9
+ ,1
+ ,32.9
+ ,33.5
+ ,32.8
+ ,1
+ ,32.9
+ ,32.9
+ ,31.9
+ ,1
+ ,32.8
+ ,32.9
+ ,30.5
+ ,1
+ ,31.9
+ ,32.8
+ ,29.2
+ ,1
+ ,30.5
+ ,31.9
+ ,28.7
+ ,1
+ ,29.2
+ ,30.5
+ ,28.4
+ ,1
+ ,28.7
+ ,29.2
+ ,28
+ ,1
+ ,28.4
+ ,28.7
+ ,27.4
+ ,1
+ ,28
+ ,28.4
+ ,26.9
+ ,1
+ ,27.4
+ ,28)
+ ,dim=c(4
+ ,59)
+ ,dimnames=list(c('Y'
+ ,'X'
+ ,'Y1'
+ ,'Y2')
+ ,1:59))
> y <- array(NA,dim=c(4,59),dimnames=list(c('Y','X','Y1','Y2'),1:59))
> 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 Y1 Y2 M1 M2 M3 M4 M5 M6 M7 M8 M9 M10 M11 t
1 54.8 0 56.0 56.6 1 0 0 0 0 0 0 0 0 0 0 1
2 52.7 0 54.8 56.0 0 1 0 0 0 0 0 0 0 0 0 2
3 50.9 0 52.7 54.8 0 0 1 0 0 0 0 0 0 0 0 3
4 50.6 0 50.9 52.7 0 0 0 1 0 0 0 0 0 0 0 4
5 52.1 0 50.6 50.9 0 0 0 0 1 0 0 0 0 0 0 5
6 53.3 0 52.1 50.6 0 0 0 0 0 1 0 0 0 0 0 6
7 53.9 0 53.3 52.1 0 0 0 0 0 0 1 0 0 0 0 7
8 54.3 0 53.9 53.3 0 0 0 0 0 0 0 1 0 0 0 8
9 54.2 0 54.3 53.9 0 0 0 0 0 0 0 0 1 0 0 9
10 54.2 0 54.2 54.3 0 0 0 0 0 0 0 0 0 1 0 10
11 53.5 0 54.2 54.2 0 0 0 0 0 0 0 0 0 0 1 11
12 51.4 0 53.5 54.2 0 0 0 0 0 0 0 0 0 0 0 12
13 50.5 0 51.4 53.5 1 0 0 0 0 0 0 0 0 0 0 13
14 50.3 0 50.5 51.4 0 1 0 0 0 0 0 0 0 0 0 14
15 49.8 0 50.3 50.5 0 0 1 0 0 0 0 0 0 0 0 15
16 50.7 0 49.8 50.3 0 0 0 1 0 0 0 0 0 0 0 16
17 52.8 0 50.7 49.8 0 0 0 0 1 0 0 0 0 0 0 17
18 55.3 0 52.8 50.7 0 0 0 0 0 1 0 0 0 0 0 18
19 57.3 0 55.3 52.8 0 0 0 0 0 0 1 0 0 0 0 19
20 57.5 0 57.3 55.3 0 0 0 0 0 0 0 1 0 0 0 20
21 56.8 0 57.5 57.3 0 0 0 0 0 0 0 0 1 0 0 21
22 56.4 0 56.8 57.5 0 0 0 0 0 0 0 0 0 1 0 22
23 56.3 0 56.4 56.8 0 0 0 0 0 0 0 0 0 0 1 23
24 56.4 0 56.3 56.4 0 0 0 0 0 0 0 0 0 0 0 24
25 57.0 0 56.4 56.3 1 0 0 0 0 0 0 0 0 0 0 25
26 57.9 0 57.0 56.4 0 1 0 0 0 0 0 0 0 0 0 26
27 58.9 0 57.9 57.0 0 0 1 0 0 0 0 0 0 0 0 27
28 58.8 0 58.9 57.9 0 0 0 1 0 0 0 0 0 0 0 28
29 56.5 1 58.8 58.9 0 0 0 0 1 0 0 0 0 0 0 29
30 51.9 1 56.5 58.8 0 0 0 0 0 1 0 0 0 0 0 30
31 47.4 1 51.9 56.5 0 0 0 0 0 0 1 0 0 0 0 31
32 44.9 1 47.4 51.9 0 0 0 0 0 0 0 1 0 0 0 32
33 43.9 1 44.9 47.4 0 0 0 0 0 0 0 0 1 0 0 33
34 43.4 1 43.9 44.9 0 0 0 0 0 0 0 0 0 1 0 34
35 42.9 1 43.4 43.9 0 0 0 0 0 0 0 0 0 0 1 35
36 42.6 1 42.9 43.4 0 0 0 0 0 0 0 0 0 0 0 36
37 42.2 1 42.6 42.9 1 0 0 0 0 0 0 0 0 0 0 37
38 41.2 1 42.2 42.6 0 1 0 0 0 0 0 0 0 0 0 38
39 40.2 1 41.2 42.2 0 0 1 0 0 0 0 0 0 0 0 39
40 39.3 1 40.2 41.2 0 0 0 1 0 0 0 0 0 0 0 40
41 38.5 1 39.3 40.2 0 0 0 0 1 0 0 0 0 0 0 41
42 38.3 1 38.5 39.3 0 0 0 0 0 1 0 0 0 0 0 42
43 37.9 1 38.3 38.5 0 0 0 0 0 0 1 0 0 0 0 43
44 37.6 1 37.9 38.3 0 0 0 0 0 0 0 1 0 0 0 44
45 37.3 1 37.6 37.9 0 0 0 0 0 0 0 0 1 0 0 45
46 36.0 1 37.3 37.6 0 0 0 0 0 0 0 0 0 1 0 46
47 34.5 1 36.0 37.3 0 0 0 0 0 0 0 0 0 0 1 47
48 33.5 1 34.5 36.0 0 0 0 0 0 0 0 0 0 0 0 48
49 32.9 1 33.5 34.5 1 0 0 0 0 0 0 0 0 0 0 49
50 32.9 1 32.9 33.5 0 1 0 0 0 0 0 0 0 0 0 50
51 32.8 1 32.9 32.9 0 0 1 0 0 0 0 0 0 0 0 51
52 31.9 1 32.8 32.9 0 0 0 1 0 0 0 0 0 0 0 52
53 30.5 1 31.9 32.8 0 0 0 0 1 0 0 0 0 0 0 53
54 29.2 1 30.5 31.9 0 0 0 0 0 1 0 0 0 0 0 54
55 28.7 1 29.2 30.5 0 0 0 0 0 0 1 0 0 0 0 55
56 28.4 1 28.7 29.2 0 0 0 0 0 0 0 1 0 0 0 56
57 28.0 1 28.4 28.7 0 0 0 0 0 0 0 0 1 0 0 57
58 27.4 1 28.0 28.4 0 0 0 0 0 0 0 0 0 1 0 58
59 26.9 1 27.4 28.0 0 0 0 0 0 0 0 0 0 0 1 59
> k <- length(x[1,])
> df <- as.data.frame(x)
> (mylm <- lm(df))
Call:
lm(formula = df)
Coefficients:
(Intercept) X Y1 Y2 M1 M2
3.460949 -1.358908 1.575945 -0.641365 0.301634 0.112538
M3 M4 M5 M6 M7 M8
0.071505 0.263299 0.460167 0.100279 0.184482 0.262352
M9 M10 M11 t
0.194356 0.104842 0.009885 -0.003196
> (mysum <- summary(mylm))
Call:
lm(formula = df)
Residuals:
Min 1Q Median 3Q Max
-1.57367 -0.49339 0.08477 0.46796 1.34984
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 3.460949 1.393426 2.484 0.0170 *
X -1.358908 0.527899 -2.574 0.0136 *
Y1 1.575945 0.106811 14.754 < 2e-16 ***
Y2 -0.641365 0.104974 -6.110 2.53e-07 ***
M1 0.301634 0.539082 0.560 0.5787
M2 0.112538 0.538915 0.209 0.8356
M3 0.071505 0.538747 0.133 0.8950
M4 0.263299 0.538875 0.489 0.6276
M5 0.460167 0.546463 0.842 0.4044
M6 0.100279 0.545736 0.184 0.8551
M7 0.184482 0.540479 0.341 0.7345
M8 0.262352 0.539049 0.487 0.6289
M9 0.194356 0.538981 0.361 0.7202
M10 0.104842 0.538855 0.195 0.8467
M11 0.009885 0.538931 0.018 0.9855
t -0.003196 0.016433 -0.194 0.8467
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 0.801 on 43 degrees of freedom
Multiple R-squared: 0.9954, Adjusted R-squared: 0.9938
F-statistic: 618.7 on 15 and 43 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.3837977 0.767595410 0.616202295
[2,] 0.5680140 0.863972021 0.431986010
[3,] 0.6815781 0.636843718 0.318421859
[4,] 0.5862693 0.827461498 0.413730749
[5,] 0.6385257 0.722948605 0.361474302
[6,] 0.9479468 0.104106380 0.052053190
[7,] 0.9268083 0.146383363 0.073191681
[8,] 0.9115660 0.176868052 0.088434026
[9,] 0.8879159 0.224168203 0.112084102
[10,] 0.9462478 0.107504348 0.053752174
[11,] 0.9968806 0.006238766 0.003119383
[12,] 0.9962503 0.007499477 0.003749739
[13,] 0.9922029 0.015594153 0.007797076
[14,] 0.9968096 0.006380777 0.003190388
[15,] 0.9956746 0.008650862 0.004325431
[16,] 0.9914409 0.017118291 0.008559146
[17,] 0.9867013 0.026597452 0.013298726
[18,] 0.9704162 0.059167696 0.029583848
[19,] 0.9465409 0.106918291 0.053459145
[20,] 0.9486931 0.102613713 0.051306856
[21,] 0.8953787 0.209242680 0.104621340
[22,] 0.7982171 0.403565737 0.201782869
> postscript(file="/var/www/html/rcomp/tmp/19nau1258649988.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/221241258649988.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/3zhq01258649988.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/4mugl1258649988.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/58ehn1258649988.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 = 59
Frequency = 1
1 2 3 4 5 6
-0.91104172 -1.31243475 -0.52835971 0.47287628 1.09753087 0.10428916
7 8 9 10 11 12
-0.30580503 -0.15640797 -0.43077501 0.07607546 -0.58990765 -1.57366588
13 14 15 16 17 18
0.08842406 0.15220036 -0.56561009 0.80549096 0.97278617 1.10361544
19 20 21 22 23 24
0.42961203 -0.99353923 -0.65480655 0.26933777 0.44891365 0.46304266
25 26 27 28 29 30
0.54287314 0.75373501 0.76443263 -0.52288208 -0.85868678 -1.53506562
31 32 33 34 35 36
-0.34186687 1.22493218 1.34984416 0.91508685 0.65984775 0.84021816
37 38 39 40 41 42
0.29388059 -0.07585877 0.28776897 0.13375050 -0.08293622 0.76367600
43 44 45 46 47 48
0.08476529 0.21219629 0.19962572 -0.72729025 -0.27281820 0.27040505
49 50 51 52 53 54
-0.01413606 0.48235815 0.04176820 -0.88923567 -1.12869404 -0.43651499
55 56 57 58 59
0.13329458 -0.28718127 -0.46388832 -0.53320982 -0.24603555
> postscript(file="/var/www/html/rcomp/tmp/6pxki1258649988.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 = 59
Frequency = 1
lag(myerror, k = 1) myerror
0 -0.91104172 NA
1 -1.31243475 -0.91104172
2 -0.52835971 -1.31243475
3 0.47287628 -0.52835971
4 1.09753087 0.47287628
5 0.10428916 1.09753087
6 -0.30580503 0.10428916
7 -0.15640797 -0.30580503
8 -0.43077501 -0.15640797
9 0.07607546 -0.43077501
10 -0.58990765 0.07607546
11 -1.57366588 -0.58990765
12 0.08842406 -1.57366588
13 0.15220036 0.08842406
14 -0.56561009 0.15220036
15 0.80549096 -0.56561009
16 0.97278617 0.80549096
17 1.10361544 0.97278617
18 0.42961203 1.10361544
19 -0.99353923 0.42961203
20 -0.65480655 -0.99353923
21 0.26933777 -0.65480655
22 0.44891365 0.26933777
23 0.46304266 0.44891365
24 0.54287314 0.46304266
25 0.75373501 0.54287314
26 0.76443263 0.75373501
27 -0.52288208 0.76443263
28 -0.85868678 -0.52288208
29 -1.53506562 -0.85868678
30 -0.34186687 -1.53506562
31 1.22493218 -0.34186687
32 1.34984416 1.22493218
33 0.91508685 1.34984416
34 0.65984775 0.91508685
35 0.84021816 0.65984775
36 0.29388059 0.84021816
37 -0.07585877 0.29388059
38 0.28776897 -0.07585877
39 0.13375050 0.28776897
40 -0.08293622 0.13375050
41 0.76367600 -0.08293622
42 0.08476529 0.76367600
43 0.21219629 0.08476529
44 0.19962572 0.21219629
45 -0.72729025 0.19962572
46 -0.27281820 -0.72729025
47 0.27040505 -0.27281820
48 -0.01413606 0.27040505
49 0.48235815 -0.01413606
50 0.04176820 0.48235815
51 -0.88923567 0.04176820
52 -1.12869404 -0.88923567
53 -0.43651499 -1.12869404
54 0.13329458 -0.43651499
55 -0.28718127 0.13329458
56 -0.46388832 -0.28718127
57 -0.53320982 -0.46388832
58 -0.24603555 -0.53320982
59 NA -0.24603555
> dum1 <- dum[2:length(myerror),]
> dum1
lag(myerror, k = 1) myerror
[1,] -1.31243475 -0.91104172
[2,] -0.52835971 -1.31243475
[3,] 0.47287628 -0.52835971
[4,] 1.09753087 0.47287628
[5,] 0.10428916 1.09753087
[6,] -0.30580503 0.10428916
[7,] -0.15640797 -0.30580503
[8,] -0.43077501 -0.15640797
[9,] 0.07607546 -0.43077501
[10,] -0.58990765 0.07607546
[11,] -1.57366588 -0.58990765
[12,] 0.08842406 -1.57366588
[13,] 0.15220036 0.08842406
[14,] -0.56561009 0.15220036
[15,] 0.80549096 -0.56561009
[16,] 0.97278617 0.80549096
[17,] 1.10361544 0.97278617
[18,] 0.42961203 1.10361544
[19,] -0.99353923 0.42961203
[20,] -0.65480655 -0.99353923
[21,] 0.26933777 -0.65480655
[22,] 0.44891365 0.26933777
[23,] 0.46304266 0.44891365
[24,] 0.54287314 0.46304266
[25,] 0.75373501 0.54287314
[26,] 0.76443263 0.75373501
[27,] -0.52288208 0.76443263
[28,] -0.85868678 -0.52288208
[29,] -1.53506562 -0.85868678
[30,] -0.34186687 -1.53506562
[31,] 1.22493218 -0.34186687
[32,] 1.34984416 1.22493218
[33,] 0.91508685 1.34984416
[34,] 0.65984775 0.91508685
[35,] 0.84021816 0.65984775
[36,] 0.29388059 0.84021816
[37,] -0.07585877 0.29388059
[38,] 0.28776897 -0.07585877
[39,] 0.13375050 0.28776897
[40,] -0.08293622 0.13375050
[41,] 0.76367600 -0.08293622
[42,] 0.08476529 0.76367600
[43,] 0.21219629 0.08476529
[44,] 0.19962572 0.21219629
[45,] -0.72729025 0.19962572
[46,] -0.27281820 -0.72729025
[47,] 0.27040505 -0.27281820
[48,] -0.01413606 0.27040505
[49,] 0.48235815 -0.01413606
[50,] 0.04176820 0.48235815
[51,] -0.88923567 0.04176820
[52,] -1.12869404 -0.88923567
[53,] -0.43651499 -1.12869404
[54,] 0.13329458 -0.43651499
[55,] -0.28718127 0.13329458
[56,] -0.46388832 -0.28718127
[57,] -0.53320982 -0.46388832
[58,] -0.24603555 -0.53320982
> z <- as.data.frame(dum1)
> z
lag(myerror, k = 1) myerror
1 -1.31243475 -0.91104172
2 -0.52835971 -1.31243475
3 0.47287628 -0.52835971
4 1.09753087 0.47287628
5 0.10428916 1.09753087
6 -0.30580503 0.10428916
7 -0.15640797 -0.30580503
8 -0.43077501 -0.15640797
9 0.07607546 -0.43077501
10 -0.58990765 0.07607546
11 -1.57366588 -0.58990765
12 0.08842406 -1.57366588
13 0.15220036 0.08842406
14 -0.56561009 0.15220036
15 0.80549096 -0.56561009
16 0.97278617 0.80549096
17 1.10361544 0.97278617
18 0.42961203 1.10361544
19 -0.99353923 0.42961203
20 -0.65480655 -0.99353923
21 0.26933777 -0.65480655
22 0.44891365 0.26933777
23 0.46304266 0.44891365
24 0.54287314 0.46304266
25 0.75373501 0.54287314
26 0.76443263 0.75373501
27 -0.52288208 0.76443263
28 -0.85868678 -0.52288208
29 -1.53506562 -0.85868678
30 -0.34186687 -1.53506562
31 1.22493218 -0.34186687
32 1.34984416 1.22493218
33 0.91508685 1.34984416
34 0.65984775 0.91508685
35 0.84021816 0.65984775
36 0.29388059 0.84021816
37 -0.07585877 0.29388059
38 0.28776897 -0.07585877
39 0.13375050 0.28776897
40 -0.08293622 0.13375050
41 0.76367600 -0.08293622
42 0.08476529 0.76367600
43 0.21219629 0.08476529
44 0.19962572 0.21219629
45 -0.72729025 0.19962572
46 -0.27281820 -0.72729025
47 0.27040505 -0.27281820
48 -0.01413606 0.27040505
49 0.48235815 -0.01413606
50 0.04176820 0.48235815
51 -0.88923567 0.04176820
52 -1.12869404 -0.88923567
53 -0.43651499 -1.12869404
54 0.13329458 -0.43651499
55 -0.28718127 0.13329458
56 -0.46388832 -0.28718127
57 -0.53320982 -0.46388832
58 -0.24603555 -0.53320982
> 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/7tbf01258649988.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/84llt1258649988.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/9yeta1258649988.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/10rarr1258649988.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/11tg7w1258649988.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/125bs71258649988.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/13s1yx1258649988.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/14p1av1258649988.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/15ymqw1258649988.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/16zh0b1258649988.tab")
+ }
>
> system("convert tmp/19nau1258649988.ps tmp/19nau1258649988.png")
> system("convert tmp/221241258649988.ps tmp/221241258649988.png")
> system("convert tmp/3zhq01258649988.ps tmp/3zhq01258649988.png")
> system("convert tmp/4mugl1258649988.ps tmp/4mugl1258649988.png")
> system("convert tmp/58ehn1258649988.ps tmp/58ehn1258649988.png")
> system("convert tmp/6pxki1258649988.ps tmp/6pxki1258649988.png")
> system("convert tmp/7tbf01258649988.ps tmp/7tbf01258649988.png")
> system("convert tmp/84llt1258649988.ps tmp/84llt1258649988.png")
> system("convert tmp/9yeta1258649988.ps tmp/9yeta1258649988.png")
> system("convert tmp/10rarr1258649988.ps tmp/10rarr1258649988.png")
>
>
> proc.time()
user system elapsed
2.408 1.592 2.803