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(37,1,30,1,47,1,35,1,30,1,43,1,82,1,40,1,47,1,19,0,52,1,136,1,80,1,42,1,54,1,66,1,81,1,63,1,137,1,72,1,107,1,58,1,36,1,52,1,79,1,77,1,54,1,84,1,48,1,96,1,83,1,66,1,61,1,53,1,30,1,74,1,69,1,59,1,42,1,65,1,70,1,100,1,63,1,105,1,82,1,81,1,75,1,102,1,121,1,98,1,76,1,77,1,63,1,37,1,35,1,23,0,40,1,29,0,37,1,51,1,20,0,28,0,13,0,22,0,25,0,13,0,16,0,13,0,16,0,17,0,9,0,17,0,25,0,14,0,8,0,7,0,10,0,7,0,10,0,3,0),dim=c(2,80),dimnames=list(c('Sol.KIA','dummy'),1:80))
> y <- array(NA,dim=c(2,80),dimnames=list(c('Sol.KIA','dummy'),1:80))
> 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
Sol.KIA dummy M1 M2 M3 M4 M5 M6 M7 M8 M9 M10 M11 t
1 37 1 1 0 0 0 0 0 0 0 0 0 0 1
2 30 1 0 1 0 0 0 0 0 0 0 0 0 2
3 47 1 0 0 1 0 0 0 0 0 0 0 0 3
4 35 1 0 0 0 1 0 0 0 0 0 0 0 4
5 30 1 0 0 0 0 1 0 0 0 0 0 0 5
6 43 1 0 0 0 0 0 1 0 0 0 0 0 6
7 82 1 0 0 0 0 0 0 1 0 0 0 0 7
8 40 1 0 0 0 0 0 0 0 1 0 0 0 8
9 47 1 0 0 0 0 0 0 0 0 1 0 0 9
10 19 0 0 0 0 0 0 0 0 0 0 1 0 10
11 52 1 0 0 0 0 0 0 0 0 0 0 1 11
12 136 1 0 0 0 0 0 0 0 0 0 0 0 12
13 80 1 1 0 0 0 0 0 0 0 0 0 0 13
14 42 1 0 1 0 0 0 0 0 0 0 0 0 14
15 54 1 0 0 1 0 0 0 0 0 0 0 0 15
16 66 1 0 0 0 1 0 0 0 0 0 0 0 16
17 81 1 0 0 0 0 1 0 0 0 0 0 0 17
18 63 1 0 0 0 0 0 1 0 0 0 0 0 18
19 137 1 0 0 0 0 0 0 1 0 0 0 0 19
20 72 1 0 0 0 0 0 0 0 1 0 0 0 20
21 107 1 0 0 0 0 0 0 0 0 1 0 0 21
22 58 1 0 0 0 0 0 0 0 0 0 1 0 22
23 36 1 0 0 0 0 0 0 0 0 0 0 1 23
24 52 1 0 0 0 0 0 0 0 0 0 0 0 24
25 79 1 1 0 0 0 0 0 0 0 0 0 0 25
26 77 1 0 1 0 0 0 0 0 0 0 0 0 26
27 54 1 0 0 1 0 0 0 0 0 0 0 0 27
28 84 1 0 0 0 1 0 0 0 0 0 0 0 28
29 48 1 0 0 0 0 1 0 0 0 0 0 0 29
30 96 1 0 0 0 0 0 1 0 0 0 0 0 30
31 83 1 0 0 0 0 0 0 1 0 0 0 0 31
32 66 1 0 0 0 0 0 0 0 1 0 0 0 32
33 61 1 0 0 0 0 0 0 0 0 1 0 0 33
34 53 1 0 0 0 0 0 0 0 0 0 1 0 34
35 30 1 0 0 0 0 0 0 0 0 0 0 1 35
36 74 1 0 0 0 0 0 0 0 0 0 0 0 36
37 69 1 1 0 0 0 0 0 0 0 0 0 0 37
38 59 1 0 1 0 0 0 0 0 0 0 0 0 38
39 42 1 0 0 1 0 0 0 0 0 0 0 0 39
40 65 1 0 0 0 1 0 0 0 0 0 0 0 40
41 70 1 0 0 0 0 1 0 0 0 0 0 0 41
42 100 1 0 0 0 0 0 1 0 0 0 0 0 42
43 63 1 0 0 0 0 0 0 1 0 0 0 0 43
44 105 1 0 0 0 0 0 0 0 1 0 0 0 44
45 82 1 0 0 0 0 0 0 0 0 1 0 0 45
46 81 1 0 0 0 0 0 0 0 0 0 1 0 46
47 75 1 0 0 0 0 0 0 0 0 0 0 1 47
48 102 1 0 0 0 0 0 0 0 0 0 0 0 48
49 121 1 1 0 0 0 0 0 0 0 0 0 0 49
50 98 1 0 1 0 0 0 0 0 0 0 0 0 50
51 76 1 0 0 1 0 0 0 0 0 0 0 0 51
52 77 1 0 0 0 1 0 0 0 0 0 0 0 52
53 63 1 0 0 0 0 1 0 0 0 0 0 0 53
54 37 1 0 0 0 0 0 1 0 0 0 0 0 54
55 35 1 0 0 0 0 0 0 1 0 0 0 0 55
56 23 0 0 0 0 0 0 0 0 1 0 0 0 56
57 40 1 0 0 0 0 0 0 0 0 1 0 0 57
58 29 0 0 0 0 0 0 0 0 0 0 1 0 58
59 37 1 0 0 0 0 0 0 0 0 0 0 1 59
60 51 1 0 0 0 0 0 0 0 0 0 0 0 60
61 20 0 1 0 0 0 0 0 0 0 0 0 0 61
62 28 0 0 1 0 0 0 0 0 0 0 0 0 62
63 13 0 0 0 1 0 0 0 0 0 0 0 0 63
64 22 0 0 0 0 1 0 0 0 0 0 0 0 64
65 25 0 0 0 0 0 1 0 0 0 0 0 0 65
66 13 0 0 0 0 0 0 1 0 0 0 0 0 66
67 16 0 0 0 0 0 0 0 1 0 0 0 0 67
68 13 0 0 0 0 0 0 0 0 1 0 0 0 68
69 16 0 0 0 0 0 0 0 0 0 1 0 0 69
70 17 0 0 0 0 0 0 0 0 0 0 1 0 70
71 9 0 0 0 0 0 0 0 0 0 0 0 1 71
72 17 0 0 0 0 0 0 0 0 0 0 0 0 72
73 25 0 1 0 0 0 0 0 0 0 0 0 0 73
74 14 0 0 1 0 0 0 0 0 0 0 0 0 74
75 8 0 0 0 1 0 0 0 0 0 0 0 0 75
76 7 0 0 0 0 1 0 0 0 0 0 0 0 76
77 10 0 0 0 0 0 1 0 0 0 0 0 0 77
78 7 0 0 0 0 0 0 1 0 0 0 0 0 78
79 10 0 0 0 0 0 0 0 1 0 0 0 0 79
80 3 0 0 0 0 0 0 0 0 1 0 0 0 80
> k <- length(x[1,])
> df <- as.data.frame(x)
> (mylm <- lm(df))
Call:
lm(formula = df)
Coefficients:
(Intercept) dummy M1 M2 M3 M4
18.4347 56.6415 -2.9279 -14.9366 -22.8024 -14.0968
M5 M6 M7 M8 M9 M10
-18.3912 -13.9712 -4.5513 -11.4684 -12.7121 -9.9831
M11 t
-32.0151 0.1515
> (mysum <- summary(mylm))
Call:
lm(formula = df)
Residuals:
Min 1Q Median 3Q Max
-43.859 -13.320 -1.639 9.378 63.596
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 18.4347 15.6602 1.177 0.2434
dummy 56.6415 8.3197 6.808 3.55e-09 ***
M1 -2.9279 12.5411 -0.233 0.8161
M2 -14.9366 12.5226 -1.193 0.2372
M3 -22.8024 12.5061 -1.823 0.0728 .
M4 -14.0968 12.4916 -1.129 0.2632
M5 -18.3912 12.4791 -1.474 0.1453
M6 -13.9712 12.4687 -1.121 0.2666
M7 -4.5513 12.4604 -0.365 0.7161
M8 -11.4684 12.5822 -0.911 0.3654
M9 -12.7121 12.9075 -0.985 0.3283
M10 -9.9831 13.2467 -0.754 0.4537
M11 -32.0151 12.8995 -2.482 0.0156 *
t 0.1515 0.1603 0.945 0.3479
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 22.34 on 66 degrees of freedom
Multiple R-squared: 0.5907, Adjusted R-squared: 0.5101
F-statistic: 7.328 on 13 and 66 DF, p-value: 1.272e-08
> 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.3819366 7.638731e-01 6.180634e-01
[2,] 0.2403731 4.807462e-01 7.596269e-01
[3,] 0.3352172 6.704344e-01 6.647828e-01
[4,] 0.2166906 4.333811e-01 7.833094e-01
[5,] 0.2861165 5.722330e-01 7.138835e-01
[6,] 0.1965770 3.931539e-01 8.034230e-01
[7,] 0.4573856 9.147713e-01 5.426144e-01
[8,] 0.9867701 2.645981e-02 1.322990e-02
[9,] 0.9781524 4.369527e-02 2.184763e-02
[10,] 0.9659368 6.812638e-02 3.406319e-02
[11,] 0.9585589 8.288223e-02 4.144111e-02
[12,] 0.9375847 1.248306e-01 6.241531e-02
[13,] 0.9476031 1.047939e-01 5.239694e-02
[14,] 0.9396318 1.207364e-01 6.036819e-02
[15,] 0.9621353 7.572942e-02 3.786471e-02
[16,] 0.9469666 1.060667e-01 5.303335e-02
[17,] 0.9431832 1.136336e-01 5.681680e-02
[18,] 0.9393191 1.213617e-01 6.068085e-02
[19,] 0.9572182 8.556362e-02 4.278181e-02
[20,] 0.9539893 9.202139e-02 4.601070e-02
[21,] 0.9561245 8.775100e-02 4.387550e-02
[22,] 0.9636272 7.274560e-02 3.637280e-02
[23,] 0.9819440 3.611194e-02 1.805597e-02
[24,] 0.9807822 3.843555e-02 1.921777e-02
[25,] 0.9760151 4.796987e-02 2.398494e-02
[26,] 0.9788017 4.239660e-02 2.119830e-02
[27,] 0.9838573 3.228542e-02 1.614271e-02
[28,] 0.9895458 2.090839e-02 1.045419e-02
[29,] 0.9849870 3.002601e-02 1.501300e-02
[30,] 0.9763878 4.722440e-02 2.361220e-02
[31,] 0.9701928 5.961430e-02 2.980715e-02
[32,] 0.9787283 4.254346e-02 2.127173e-02
[33,] 0.9988870 2.225904e-03 1.112952e-03
[34,] 0.9998586 2.827436e-04 1.413718e-04
[35,] 0.9999716 5.678599e-05 2.839300e-05
[36,] 0.9999995 9.575477e-07 4.787739e-07
[37,] 0.9999999 2.414881e-07 1.207441e-07
[38,] 0.9999997 5.672296e-07 2.836148e-07
[39,] 0.9999998 3.218044e-07 1.609022e-07
[40,] 0.9999991 1.816827e-06 9.084137e-07
[41,] 0.9999972 5.653977e-06 2.826989e-06
[42,] 0.9999852 2.951268e-05 1.475634e-05
[43,] 0.9999326 1.348566e-04 6.742831e-05
[44,] 0.9996907 6.185627e-04 3.092814e-04
[45,] 0.9998501 2.997125e-04 1.498563e-04
[46,] 0.9991457 1.708535e-03 8.542674e-04
[47,] 0.9956578 8.684488e-03 4.342244e-03
> postscript(file="/var/www/html/rcomp/tmp/15gir1291029063.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/25gir1291029063.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/3g70u1291029063.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/4g70u1291029063.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/5g70u1291029063.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 = 80
Frequency = 1
1 2 3 4 5 6
-35.2998652 -30.4427224 -5.7284367 -26.5855795 -27.4427224 -19.0141509
7 8 9 10 11 12
10.4144205 -24.8200809 -16.7278751 9.0331312 7.2721249 59.1054582
13 14 15 16 17 18
5.8818509 -20.2610063 -0.5467206 2.5961366 21.7389937 -0.8324349
19 20 21 22 23 24
63.5961366 5.3616352 41.4538410 -10.4266622 -10.5461590 -26.7128257
25 26 27 28 29 30
3.0635669 12.9207098 -2.3650045 18.7778527 -13.0792902 30.3492812
31 32 33 34 35 36
7.7778527 -2.4566487 -6.3644429 -17.2449461 -18.3644429 -6.5311096
37 38 39 40 41 42
-8.7547170 -6.8975741 -16.1832884 -2.0404313 7.1024259 32.5309973
43 44 45 46 47 48
-14.0404313 34.7250674 12.8172731 8.9367700 24.8172731 19.6506065
49 50 51 52 53 54
41.4269991 30.2841420 15.9984277 8.1412848 -1.7158580 -32.2872866
55 56 57 58 59 60
-43.8587152 7.5482929 -31.0010108 11.7599955 -15.0010108 -33.1676774
61 62 63 64 65 66
-4.7497754 15.1073675 7.8216532 7.9645103 15.1073675 -1.4640611
67 68 69 70 71 72
-8.0354897 -4.2699910 -0.1777853 -2.0582884 11.8222147 -12.3444519
73 74 75 76 77 78
-1.5680593 -0.7109164 1.0033693 -8.8537736 -1.7109164 -9.2823450
79 80
-15.8537736 -16.0882749
> postscript(file="/var/www/html/rcomp/tmp/69gzf1291029063.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 = 80
Frequency = 1
lag(myerror, k = 1) myerror
0 -35.2998652 NA
1 -30.4427224 -35.2998652
2 -5.7284367 -30.4427224
3 -26.5855795 -5.7284367
4 -27.4427224 -26.5855795
5 -19.0141509 -27.4427224
6 10.4144205 -19.0141509
7 -24.8200809 10.4144205
8 -16.7278751 -24.8200809
9 9.0331312 -16.7278751
10 7.2721249 9.0331312
11 59.1054582 7.2721249
12 5.8818509 59.1054582
13 -20.2610063 5.8818509
14 -0.5467206 -20.2610063
15 2.5961366 -0.5467206
16 21.7389937 2.5961366
17 -0.8324349 21.7389937
18 63.5961366 -0.8324349
19 5.3616352 63.5961366
20 41.4538410 5.3616352
21 -10.4266622 41.4538410
22 -10.5461590 -10.4266622
23 -26.7128257 -10.5461590
24 3.0635669 -26.7128257
25 12.9207098 3.0635669
26 -2.3650045 12.9207098
27 18.7778527 -2.3650045
28 -13.0792902 18.7778527
29 30.3492812 -13.0792902
30 7.7778527 30.3492812
31 -2.4566487 7.7778527
32 -6.3644429 -2.4566487
33 -17.2449461 -6.3644429
34 -18.3644429 -17.2449461
35 -6.5311096 -18.3644429
36 -8.7547170 -6.5311096
37 -6.8975741 -8.7547170
38 -16.1832884 -6.8975741
39 -2.0404313 -16.1832884
40 7.1024259 -2.0404313
41 32.5309973 7.1024259
42 -14.0404313 32.5309973
43 34.7250674 -14.0404313
44 12.8172731 34.7250674
45 8.9367700 12.8172731
46 24.8172731 8.9367700
47 19.6506065 24.8172731
48 41.4269991 19.6506065
49 30.2841420 41.4269991
50 15.9984277 30.2841420
51 8.1412848 15.9984277
52 -1.7158580 8.1412848
53 -32.2872866 -1.7158580
54 -43.8587152 -32.2872866
55 7.5482929 -43.8587152
56 -31.0010108 7.5482929
57 11.7599955 -31.0010108
58 -15.0010108 11.7599955
59 -33.1676774 -15.0010108
60 -4.7497754 -33.1676774
61 15.1073675 -4.7497754
62 7.8216532 15.1073675
63 7.9645103 7.8216532
64 15.1073675 7.9645103
65 -1.4640611 15.1073675
66 -8.0354897 -1.4640611
67 -4.2699910 -8.0354897
68 -0.1777853 -4.2699910
69 -2.0582884 -0.1777853
70 11.8222147 -2.0582884
71 -12.3444519 11.8222147
72 -1.5680593 -12.3444519
73 -0.7109164 -1.5680593
74 1.0033693 -0.7109164
75 -8.8537736 1.0033693
76 -1.7109164 -8.8537736
77 -9.2823450 -1.7109164
78 -15.8537736 -9.2823450
79 -16.0882749 -15.8537736
80 NA -16.0882749
> dum1 <- dum[2:length(myerror),]
> dum1
lag(myerror, k = 1) myerror
[1,] -30.4427224 -35.2998652
[2,] -5.7284367 -30.4427224
[3,] -26.5855795 -5.7284367
[4,] -27.4427224 -26.5855795
[5,] -19.0141509 -27.4427224
[6,] 10.4144205 -19.0141509
[7,] -24.8200809 10.4144205
[8,] -16.7278751 -24.8200809
[9,] 9.0331312 -16.7278751
[10,] 7.2721249 9.0331312
[11,] 59.1054582 7.2721249
[12,] 5.8818509 59.1054582
[13,] -20.2610063 5.8818509
[14,] -0.5467206 -20.2610063
[15,] 2.5961366 -0.5467206
[16,] 21.7389937 2.5961366
[17,] -0.8324349 21.7389937
[18,] 63.5961366 -0.8324349
[19,] 5.3616352 63.5961366
[20,] 41.4538410 5.3616352
[21,] -10.4266622 41.4538410
[22,] -10.5461590 -10.4266622
[23,] -26.7128257 -10.5461590
[24,] 3.0635669 -26.7128257
[25,] 12.9207098 3.0635669
[26,] -2.3650045 12.9207098
[27,] 18.7778527 -2.3650045
[28,] -13.0792902 18.7778527
[29,] 30.3492812 -13.0792902
[30,] 7.7778527 30.3492812
[31,] -2.4566487 7.7778527
[32,] -6.3644429 -2.4566487
[33,] -17.2449461 -6.3644429
[34,] -18.3644429 -17.2449461
[35,] -6.5311096 -18.3644429
[36,] -8.7547170 -6.5311096
[37,] -6.8975741 -8.7547170
[38,] -16.1832884 -6.8975741
[39,] -2.0404313 -16.1832884
[40,] 7.1024259 -2.0404313
[41,] 32.5309973 7.1024259
[42,] -14.0404313 32.5309973
[43,] 34.7250674 -14.0404313
[44,] 12.8172731 34.7250674
[45,] 8.9367700 12.8172731
[46,] 24.8172731 8.9367700
[47,] 19.6506065 24.8172731
[48,] 41.4269991 19.6506065
[49,] 30.2841420 41.4269991
[50,] 15.9984277 30.2841420
[51,] 8.1412848 15.9984277
[52,] -1.7158580 8.1412848
[53,] -32.2872866 -1.7158580
[54,] -43.8587152 -32.2872866
[55,] 7.5482929 -43.8587152
[56,] -31.0010108 7.5482929
[57,] 11.7599955 -31.0010108
[58,] -15.0010108 11.7599955
[59,] -33.1676774 -15.0010108
[60,] -4.7497754 -33.1676774
[61,] 15.1073675 -4.7497754
[62,] 7.8216532 15.1073675
[63,] 7.9645103 7.8216532
[64,] 15.1073675 7.9645103
[65,] -1.4640611 15.1073675
[66,] -8.0354897 -1.4640611
[67,] -4.2699910 -8.0354897
[68,] -0.1777853 -4.2699910
[69,] -2.0582884 -0.1777853
[70,] 11.8222147 -2.0582884
[71,] -12.3444519 11.8222147
[72,] -1.5680593 -12.3444519
[73,] -0.7109164 -1.5680593
[74,] 1.0033693 -0.7109164
[75,] -8.8537736 1.0033693
[76,] -1.7109164 -8.8537736
[77,] -9.2823450 -1.7109164
[78,] -15.8537736 -9.2823450
[79,] -16.0882749 -15.8537736
> z <- as.data.frame(dum1)
> z
lag(myerror, k = 1) myerror
1 -30.4427224 -35.2998652
2 -5.7284367 -30.4427224
3 -26.5855795 -5.7284367
4 -27.4427224 -26.5855795
5 -19.0141509 -27.4427224
6 10.4144205 -19.0141509
7 -24.8200809 10.4144205
8 -16.7278751 -24.8200809
9 9.0331312 -16.7278751
10 7.2721249 9.0331312
11 59.1054582 7.2721249
12 5.8818509 59.1054582
13 -20.2610063 5.8818509
14 -0.5467206 -20.2610063
15 2.5961366 -0.5467206
16 21.7389937 2.5961366
17 -0.8324349 21.7389937
18 63.5961366 -0.8324349
19 5.3616352 63.5961366
20 41.4538410 5.3616352
21 -10.4266622 41.4538410
22 -10.5461590 -10.4266622
23 -26.7128257 -10.5461590
24 3.0635669 -26.7128257
25 12.9207098 3.0635669
26 -2.3650045 12.9207098
27 18.7778527 -2.3650045
28 -13.0792902 18.7778527
29 30.3492812 -13.0792902
30 7.7778527 30.3492812
31 -2.4566487 7.7778527
32 -6.3644429 -2.4566487
33 -17.2449461 -6.3644429
34 -18.3644429 -17.2449461
35 -6.5311096 -18.3644429
36 -8.7547170 -6.5311096
37 -6.8975741 -8.7547170
38 -16.1832884 -6.8975741
39 -2.0404313 -16.1832884
40 7.1024259 -2.0404313
41 32.5309973 7.1024259
42 -14.0404313 32.5309973
43 34.7250674 -14.0404313
44 12.8172731 34.7250674
45 8.9367700 12.8172731
46 24.8172731 8.9367700
47 19.6506065 24.8172731
48 41.4269991 19.6506065
49 30.2841420 41.4269991
50 15.9984277 30.2841420
51 8.1412848 15.9984277
52 -1.7158580 8.1412848
53 -32.2872866 -1.7158580
54 -43.8587152 -32.2872866
55 7.5482929 -43.8587152
56 -31.0010108 7.5482929
57 11.7599955 -31.0010108
58 -15.0010108 11.7599955
59 -33.1676774 -15.0010108
60 -4.7497754 -33.1676774
61 15.1073675 -4.7497754
62 7.8216532 15.1073675
63 7.9645103 7.8216532
64 15.1073675 7.9645103
65 -1.4640611 15.1073675
66 -8.0354897 -1.4640611
67 -4.2699910 -8.0354897
68 -0.1777853 -4.2699910
69 -2.0582884 -0.1777853
70 11.8222147 -2.0582884
71 -12.3444519 11.8222147
72 -1.5680593 -12.3444519
73 -0.7109164 -1.5680593
74 1.0033693 -0.7109164
75 -8.8537736 1.0033693
76 -1.7109164 -8.8537736
77 -9.2823450 -1.7109164
78 -15.8537736 -9.2823450
79 -16.0882749 -15.8537736
> 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/71qyi1291029063.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/81qyi1291029063.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/91qyi1291029063.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/10czy31291029063.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/11xher1291029063.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/1210dx1291029063.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/13fsan1291029063.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/140s9b1291029063.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/154bpz1291029063.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/167t651291029063.tab")
+ }
>
> try(system("convert tmp/15gir1291029063.ps tmp/15gir1291029063.png",intern=TRUE))
character(0)
> try(system("convert tmp/25gir1291029063.ps tmp/25gir1291029063.png",intern=TRUE))
character(0)
> try(system("convert tmp/3g70u1291029063.ps tmp/3g70u1291029063.png",intern=TRUE))
character(0)
> try(system("convert tmp/4g70u1291029063.ps tmp/4g70u1291029063.png",intern=TRUE))
character(0)
> try(system("convert tmp/5g70u1291029063.ps tmp/5g70u1291029063.png",intern=TRUE))
character(0)
> try(system("convert tmp/69gzf1291029063.ps tmp/69gzf1291029063.png",intern=TRUE))
character(0)
> try(system("convert tmp/71qyi1291029063.ps tmp/71qyi1291029063.png",intern=TRUE))
character(0)
> try(system("convert tmp/81qyi1291029063.ps tmp/81qyi1291029063.png",intern=TRUE))
character(0)
> try(system("convert tmp/91qyi1291029063.ps tmp/91qyi1291029063.png",intern=TRUE))
character(0)
> try(system("convert tmp/10czy31291029063.ps tmp/10czy31291029063.png",intern=TRUE))
character(0)
>
>
> proc.time()
user system elapsed
2.666 1.615 6.429