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(105.5073942
+ ,1
+ ,95.84395716
+ ,100
+ ,118.1540031
+ ,1
+ ,105.5073942
+ ,95.84395716
+ ,101.8612953
+ ,1
+ ,118.1540031
+ ,105.5073942
+ ,109.8419174
+ ,1
+ ,101.8612953
+ ,118.1540031
+ ,105.6348802
+ ,1
+ ,109.8419174
+ ,101.8612953
+ ,112.927078
+ ,1
+ ,105.6348802
+ ,109.8419174
+ ,133.0698623
+ ,1
+ ,112.927078
+ ,105.6348802
+ ,125.6756757
+ ,1
+ ,133.0698623
+ ,112.927078
+ ,146.736359
+ ,1
+ ,125.6756757
+ ,133.0698623
+ ,142.5803162
+ ,1
+ ,146.736359
+ ,125.6756757
+ ,106.1448241
+ ,1
+ ,142.5803162
+ ,146.736359
+ ,126.5170831
+ ,1
+ ,106.1448241
+ ,142.5803162
+ ,132.7893932
+ ,1
+ ,126.5170831
+ ,106.1448241
+ ,121.2391637
+ ,1
+ ,132.7893932
+ ,126.5170831
+ ,114.5079041
+ ,1
+ ,121.2391637
+ ,132.7893932
+ ,146.1499235
+ ,1
+ ,114.5079041
+ ,121.2391637
+ ,146.1244263
+ ,1
+ ,146.1499235
+ ,114.5079041
+ ,128.5058644
+ ,1
+ ,146.1244263
+ ,146.1499235
+ ,155.5838858
+ ,1
+ ,128.5058644
+ ,146.1244263
+ ,125.0382458
+ ,1
+ ,155.5838858
+ ,128.5058644
+ ,136.8944416
+ ,1
+ ,125.0382458
+ ,155.5838858
+ ,142.2233554
+ ,1
+ ,136.8944416
+ ,125.0382458
+ ,117.7715451
+ ,1
+ ,142.2233554
+ ,136.8944416
+ ,120.627231
+ ,1
+ ,117.7715451
+ ,142.2233554
+ ,127.7664457
+ ,1
+ ,120.627231
+ ,117.7715451
+ ,135.1096379
+ ,1
+ ,127.7664457
+ ,120.627231
+ ,105.7113717
+ ,1
+ ,135.1096379
+ ,127.7664457
+ ,117.9245283
+ ,1
+ ,105.7113717
+ ,135.1096379
+ ,120.754717
+ ,1
+ ,117.9245283
+ ,105.7113717
+ ,107.572667
+ ,1
+ ,120.754717
+ ,117.9245283
+ ,130.4436512
+ ,1
+ ,107.572667
+ ,120.754717
+ ,107.2157063
+ ,1
+ ,130.4436512
+ ,107.572667
+ ,105.0739419
+ ,1
+ ,107.2157063
+ ,130.4436512
+ ,130.1121877
+ ,1
+ ,105.0739419
+ ,107.2157063
+ ,109.6379398
+ ,1
+ ,130.1121877
+ ,105.0739419
+ ,116.7261601
+ ,1
+ ,109.6379398
+ ,130.1121877
+ ,97.11881693
+ ,0
+ ,116.7261601
+ ,109.6379398
+ ,140.8975013
+ ,1
+ ,97.11881693
+ ,116.7261601
+ ,108.2865885
+ ,1
+ ,140.8975013
+ ,97.11881693
+ ,97.65425803
+ ,0
+ ,108.2865885
+ ,140.8975013
+ ,112.0346762
+ ,1
+ ,97.65425803
+ ,108.2865885
+ ,123.0494646
+ ,1
+ ,112.0346762
+ ,97.65425803
+ ,112.4171341
+ ,1
+ ,123.0494646
+ ,112.0346762
+ ,116.4966854
+ ,1
+ ,112.4171341
+ ,123.0494646
+ ,104.6914839
+ ,1
+ ,116.4966854
+ ,112.4171341
+ ,122.2335543
+ ,1
+ ,104.6914839
+ ,116.4966854
+ ,99.79602244
+ ,0
+ ,122.2335543
+ ,104.6914839
+ ,96.71086181
+ ,0
+ ,99.79602244
+ ,122.2335543
+ ,112.3151453
+ ,1
+ ,96.71086181
+ ,99.79602244
+ ,102.5497195
+ ,1
+ ,112.3151453
+ ,96.71086181
+ ,104.5385008
+ ,1
+ ,102.5497195
+ ,112.3151453
+ ,122.0805711
+ ,1
+ ,104.5385008
+ ,102.5497195
+ ,80.64762876
+ ,0
+ ,122.0805711
+ ,104.5385008
+ ,91.40744518
+ ,0
+ ,80.64762876
+ ,122.0805711
+ ,99.51555329
+ ,0
+ ,91.40744518
+ ,80.64762876
+ ,106.527282
+ ,1
+ ,99.51555329
+ ,91.40744518
+ ,98.49566548
+ ,0
+ ,106.527282
+ ,99.51555329
+ ,106.7567568
+ ,1
+ ,98.49566548
+ ,106.527282)
+ ,dim=c(4
+ ,58)
+ ,dimnames=list(c('Y'
+ ,'X'
+ ,'Y1'
+ ,'Y2')
+ ,1:58))
> y <- array(NA,dim=c(4,58),dimnames=list(c('Y','X','Y1','Y2'),1:58))
> 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 105.50739 1 95.84396 100.00000 1 0 0 0 0 0 0 0 0 0 0 1
2 118.15400 1 105.50739 95.84396 0 1 0 0 0 0 0 0 0 0 0 2
3 101.86130 1 118.15400 105.50739 0 0 1 0 0 0 0 0 0 0 0 3
4 109.84192 1 101.86130 118.15400 0 0 0 1 0 0 0 0 0 0 0 4
5 105.63488 1 109.84192 101.86130 0 0 0 0 1 0 0 0 0 0 0 5
6 112.92708 1 105.63488 109.84192 0 0 0 0 0 1 0 0 0 0 0 6
7 133.06986 1 112.92708 105.63488 0 0 0 0 0 0 1 0 0 0 0 7
8 125.67568 1 133.06986 112.92708 0 0 0 0 0 0 0 1 0 0 0 8
9 146.73636 1 125.67568 133.06986 0 0 0 0 0 0 0 0 1 0 0 9
10 142.58032 1 146.73636 125.67568 0 0 0 0 0 0 0 0 0 1 0 10
11 106.14482 1 142.58032 146.73636 0 0 0 0 0 0 0 0 0 0 1 11
12 126.51708 1 106.14482 142.58032 0 0 0 0 0 0 0 0 0 0 0 12
13 132.78939 1 126.51708 106.14482 1 0 0 0 0 0 0 0 0 0 0 13
14 121.23916 1 132.78939 126.51708 0 1 0 0 0 0 0 0 0 0 0 14
15 114.50790 1 121.23916 132.78939 0 0 1 0 0 0 0 0 0 0 0 15
16 146.14992 1 114.50790 121.23916 0 0 0 1 0 0 0 0 0 0 0 16
17 146.12443 1 146.14992 114.50790 0 0 0 0 1 0 0 0 0 0 0 17
18 128.50586 1 146.12443 146.14992 0 0 0 0 0 1 0 0 0 0 0 18
19 155.58389 1 128.50586 146.12443 0 0 0 0 0 0 1 0 0 0 0 19
20 125.03825 1 155.58389 128.50586 0 0 0 0 0 0 0 1 0 0 0 20
21 136.89444 1 125.03825 155.58389 0 0 0 0 0 0 0 0 1 0 0 21
22 142.22336 1 136.89444 125.03825 0 0 0 0 0 0 0 0 0 1 0 22
23 117.77155 1 142.22336 136.89444 0 0 0 0 0 0 0 0 0 0 1 23
24 120.62723 1 117.77155 142.22336 0 0 0 0 0 0 0 0 0 0 0 24
25 127.76645 1 120.62723 117.77155 1 0 0 0 0 0 0 0 0 0 0 25
26 135.10964 1 127.76645 120.62723 0 1 0 0 0 0 0 0 0 0 0 26
27 105.71137 1 135.10964 127.76645 0 0 1 0 0 0 0 0 0 0 0 27
28 117.92453 1 105.71137 135.10964 0 0 0 1 0 0 0 0 0 0 0 28
29 120.75472 1 117.92453 105.71137 0 0 0 0 1 0 0 0 0 0 0 29
30 107.57267 1 120.75472 117.92453 0 0 0 0 0 1 0 0 0 0 0 30
31 130.44365 1 107.57267 120.75472 0 0 0 0 0 0 1 0 0 0 0 31
32 107.21571 1 130.44365 107.57267 0 0 0 0 0 0 0 1 0 0 0 32
33 105.07394 1 107.21571 130.44365 0 0 0 0 0 0 0 0 1 0 0 33
34 130.11219 1 105.07394 107.21571 0 0 0 0 0 0 0 0 0 1 0 34
35 109.63794 1 130.11219 105.07394 0 0 0 0 0 0 0 0 0 0 1 35
36 116.72616 1 109.63794 130.11219 0 0 0 0 0 0 0 0 0 0 0 36
37 97.11882 0 116.72616 109.63794 1 0 0 0 0 0 0 0 0 0 0 37
38 140.89750 1 97.11882 116.72616 0 1 0 0 0 0 0 0 0 0 0 38
39 108.28659 1 140.89750 97.11882 0 0 1 0 0 0 0 0 0 0 0 39
40 97.65426 0 108.28659 140.89750 0 0 0 1 0 0 0 0 0 0 0 40
41 112.03468 1 97.65426 108.28659 0 0 0 0 1 0 0 0 0 0 0 41
42 123.04946 1 112.03468 97.65426 0 0 0 0 0 1 0 0 0 0 0 42
43 112.41713 1 123.04946 112.03468 0 0 0 0 0 0 1 0 0 0 0 43
44 116.49669 1 112.41713 123.04946 0 0 0 0 0 0 0 1 0 0 0 44
45 104.69148 1 116.49669 112.41713 0 0 0 0 0 0 0 0 1 0 0 45
46 122.23355 1 104.69148 116.49669 0 0 0 0 0 0 0 0 0 1 0 46
47 99.79602 0 122.23355 104.69148 0 0 0 0 0 0 0 0 0 0 1 47
48 96.71086 0 99.79602 122.23355 0 0 0 0 0 0 0 0 0 0 0 48
49 112.31515 1 96.71086 99.79602 1 0 0 0 0 0 0 0 0 0 0 49
50 102.54972 1 112.31515 96.71086 0 1 0 0 0 0 0 0 0 0 0 50
51 104.53850 1 102.54972 112.31515 0 0 1 0 0 0 0 0 0 0 0 51
52 122.08057 1 104.53850 102.54972 0 0 0 1 0 0 0 0 0 0 0 52
53 80.64763 0 122.08057 104.53850 0 0 0 0 1 0 0 0 0 0 0 53
54 91.40745 0 80.64763 122.08057 0 0 0 0 0 1 0 0 0 0 0 54
55 99.51555 0 91.40745 80.64763 0 0 0 0 0 0 1 0 0 0 0 55
56 106.52728 1 99.51555 91.40745 0 0 0 0 0 0 0 1 0 0 0 56
57 98.49567 0 106.52728 99.51555 0 0 0 0 0 0 0 0 1 0 0 57
58 106.75676 1 98.49567 106.52728 0 0 0 0 0 0 0 0 0 1 0 58
> k <- length(x[1,])
> df <- as.data.frame(x)
> (mylm <- lm(df))
Call:
lm(formula = df)
Coefficients:
(Intercept) X Y1 Y2 M1 M2
51.8114 18.2085 0.2712 0.1787 2.5627 5.6758
M3 M4 M5 M6 M7 M8
-13.7940 4.7073 -1.0779 -1.8577 12.8907 -4.2458
M9 M10 M11 t
2.0080 10.0896 -12.0141 -0.1236
> (mysum <- summary(mylm))
Call:
lm(formula = df)
Residuals:
Min 1Q Median 3Q Max
-18.577 -6.534 1.481 4.935 22.695
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 51.8114 20.1630 2.570 0.013820 *
X 18.2085 5.1382 3.544 0.000983 ***
Y1 0.2712 0.1286 2.108 0.041019 *
Y2 0.1787 0.1276 1.401 0.168479
M1 2.5627 8.4231 0.304 0.762439
M2 5.6758 8.3476 0.680 0.500281
M3 -13.7940 8.4555 -1.631 0.110289
M4 4.7073 7.6338 0.617 0.540803
M5 -1.0779 8.6109 -0.125 0.900978
M6 -1.8577 7.8544 -0.237 0.814181
M7 12.8907 8.0792 1.596 0.118089
M8 -4.2458 8.7416 -0.486 0.629704
M9 2.0080 7.7179 0.260 0.796000
M10 10.0896 8.2556 1.222 0.228466
M11 -12.0141 8.9022 -1.350 0.184387
t -0.1236 0.1139 -1.085 0.284301
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 11.19 on 42 degrees of freedom
Multiple R-squared: 0.6324, Adjusted R-squared: 0.5011
F-statistic: 4.816 on 15 and 42 DF, p-value: 2.735e-05
> 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.6300202 0.7399596 0.3699798
[2,] 0.8209250 0.3581501 0.1790750
[3,] 0.8949869 0.2100262 0.1050131
[4,] 0.8853892 0.2292217 0.1146108
[5,] 0.8162879 0.3674242 0.1837121
[6,] 0.8443285 0.3113430 0.1556715
[7,] 0.7878200 0.4243600 0.2121800
[8,] 0.7325944 0.5348113 0.2674056
[9,] 0.7456289 0.5087422 0.2543711
[10,] 0.7213634 0.5572733 0.2786366
[11,] 0.6309559 0.7380883 0.3690441
[12,] 0.6182693 0.7634614 0.3817307
[13,] 0.5429378 0.9141245 0.4570622
[14,] 0.5096743 0.9806514 0.4903257
[15,] 0.6343800 0.7312401 0.3656200
[16,] 0.5984901 0.8030199 0.4015099
[17,] 0.6001133 0.7997733 0.3998867
[18,] 0.4796999 0.9593999 0.5203001
[19,] 0.3601281 0.7202561 0.6398719
[20,] 0.7385027 0.5229947 0.2614973
[21,] 0.6065960 0.7868081 0.3934040
> postscript(file="/var/www/html/rcomp/tmp/1sa701258720956.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/2f8s81258720956.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/36scx1258720956.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/47qdz1258720956.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/5i64d1258720956.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 = 58
Frequency = 1
1 2 3 4 5 6
-10.8197886 -3.0406132 -4.8970833 -13.1358834 -10.6864047 -2.7763184
7 8 9 10 11 12
1.5159424 4.6154301 17.9509435 1.4466828 -15.3988634 3.7074765
13 14 15 16 17 18
8.5280146 -11.3541740 3.5193844 20.6738031 19.1785638 -3.1854720
19 20 21 22 23 24
14.0506763 -3.4296035 5.7407289 5.3559164 -0.4331428 -3.7888392
25 26 27 28 29 30
4.5073009 6.4143660 -6.6581342 -6.1620881 4.5192202 -10.7100442
31 32 33 34 35 36
0.6054061 -9.2091753 -15.2694671 6.5434714 1.8886133 -1.8361904
37 38 39 40 41 42
-3.9369338 22.6945068 1.3083924 -8.4737534 2.3194210 12.2378944
43 44 45 46 47 48
-18.5769239 3.6774620 -13.4639077 -1.4073105 13.9433930 1.9175530
49 50 51 52 53 54
1.7214069 -14.7140856 6.7274408 7.0979219 -15.3308004 4.4339402
55 56 57 58
2.4048990 4.3458867 5.0417024 -11.9387601
> postscript(file="/var/www/html/rcomp/tmp/6g06e1258720956.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 = 58
Frequency = 1
lag(myerror, k = 1) myerror
0 -10.8197886 NA
1 -3.0406132 -10.8197886
2 -4.8970833 -3.0406132
3 -13.1358834 -4.8970833
4 -10.6864047 -13.1358834
5 -2.7763184 -10.6864047
6 1.5159424 -2.7763184
7 4.6154301 1.5159424
8 17.9509435 4.6154301
9 1.4466828 17.9509435
10 -15.3988634 1.4466828
11 3.7074765 -15.3988634
12 8.5280146 3.7074765
13 -11.3541740 8.5280146
14 3.5193844 -11.3541740
15 20.6738031 3.5193844
16 19.1785638 20.6738031
17 -3.1854720 19.1785638
18 14.0506763 -3.1854720
19 -3.4296035 14.0506763
20 5.7407289 -3.4296035
21 5.3559164 5.7407289
22 -0.4331428 5.3559164
23 -3.7888392 -0.4331428
24 4.5073009 -3.7888392
25 6.4143660 4.5073009
26 -6.6581342 6.4143660
27 -6.1620881 -6.6581342
28 4.5192202 -6.1620881
29 -10.7100442 4.5192202
30 0.6054061 -10.7100442
31 -9.2091753 0.6054061
32 -15.2694671 -9.2091753
33 6.5434714 -15.2694671
34 1.8886133 6.5434714
35 -1.8361904 1.8886133
36 -3.9369338 -1.8361904
37 22.6945068 -3.9369338
38 1.3083924 22.6945068
39 -8.4737534 1.3083924
40 2.3194210 -8.4737534
41 12.2378944 2.3194210
42 -18.5769239 12.2378944
43 3.6774620 -18.5769239
44 -13.4639077 3.6774620
45 -1.4073105 -13.4639077
46 13.9433930 -1.4073105
47 1.9175530 13.9433930
48 1.7214069 1.9175530
49 -14.7140856 1.7214069
50 6.7274408 -14.7140856
51 7.0979219 6.7274408
52 -15.3308004 7.0979219
53 4.4339402 -15.3308004
54 2.4048990 4.4339402
55 4.3458867 2.4048990
56 5.0417024 4.3458867
57 -11.9387601 5.0417024
58 NA -11.9387601
> dum1 <- dum[2:length(myerror),]
> dum1
lag(myerror, k = 1) myerror
[1,] -3.0406132 -10.8197886
[2,] -4.8970833 -3.0406132
[3,] -13.1358834 -4.8970833
[4,] -10.6864047 -13.1358834
[5,] -2.7763184 -10.6864047
[6,] 1.5159424 -2.7763184
[7,] 4.6154301 1.5159424
[8,] 17.9509435 4.6154301
[9,] 1.4466828 17.9509435
[10,] -15.3988634 1.4466828
[11,] 3.7074765 -15.3988634
[12,] 8.5280146 3.7074765
[13,] -11.3541740 8.5280146
[14,] 3.5193844 -11.3541740
[15,] 20.6738031 3.5193844
[16,] 19.1785638 20.6738031
[17,] -3.1854720 19.1785638
[18,] 14.0506763 -3.1854720
[19,] -3.4296035 14.0506763
[20,] 5.7407289 -3.4296035
[21,] 5.3559164 5.7407289
[22,] -0.4331428 5.3559164
[23,] -3.7888392 -0.4331428
[24,] 4.5073009 -3.7888392
[25,] 6.4143660 4.5073009
[26,] -6.6581342 6.4143660
[27,] -6.1620881 -6.6581342
[28,] 4.5192202 -6.1620881
[29,] -10.7100442 4.5192202
[30,] 0.6054061 -10.7100442
[31,] -9.2091753 0.6054061
[32,] -15.2694671 -9.2091753
[33,] 6.5434714 -15.2694671
[34,] 1.8886133 6.5434714
[35,] -1.8361904 1.8886133
[36,] -3.9369338 -1.8361904
[37,] 22.6945068 -3.9369338
[38,] 1.3083924 22.6945068
[39,] -8.4737534 1.3083924
[40,] 2.3194210 -8.4737534
[41,] 12.2378944 2.3194210
[42,] -18.5769239 12.2378944
[43,] 3.6774620 -18.5769239
[44,] -13.4639077 3.6774620
[45,] -1.4073105 -13.4639077
[46,] 13.9433930 -1.4073105
[47,] 1.9175530 13.9433930
[48,] 1.7214069 1.9175530
[49,] -14.7140856 1.7214069
[50,] 6.7274408 -14.7140856
[51,] 7.0979219 6.7274408
[52,] -15.3308004 7.0979219
[53,] 4.4339402 -15.3308004
[54,] 2.4048990 4.4339402
[55,] 4.3458867 2.4048990
[56,] 5.0417024 4.3458867
[57,] -11.9387601 5.0417024
> z <- as.data.frame(dum1)
> z
lag(myerror, k = 1) myerror
1 -3.0406132 -10.8197886
2 -4.8970833 -3.0406132
3 -13.1358834 -4.8970833
4 -10.6864047 -13.1358834
5 -2.7763184 -10.6864047
6 1.5159424 -2.7763184
7 4.6154301 1.5159424
8 17.9509435 4.6154301
9 1.4466828 17.9509435
10 -15.3988634 1.4466828
11 3.7074765 -15.3988634
12 8.5280146 3.7074765
13 -11.3541740 8.5280146
14 3.5193844 -11.3541740
15 20.6738031 3.5193844
16 19.1785638 20.6738031
17 -3.1854720 19.1785638
18 14.0506763 -3.1854720
19 -3.4296035 14.0506763
20 5.7407289 -3.4296035
21 5.3559164 5.7407289
22 -0.4331428 5.3559164
23 -3.7888392 -0.4331428
24 4.5073009 -3.7888392
25 6.4143660 4.5073009
26 -6.6581342 6.4143660
27 -6.1620881 -6.6581342
28 4.5192202 -6.1620881
29 -10.7100442 4.5192202
30 0.6054061 -10.7100442
31 -9.2091753 0.6054061
32 -15.2694671 -9.2091753
33 6.5434714 -15.2694671
34 1.8886133 6.5434714
35 -1.8361904 1.8886133
36 -3.9369338 -1.8361904
37 22.6945068 -3.9369338
38 1.3083924 22.6945068
39 -8.4737534 1.3083924
40 2.3194210 -8.4737534
41 12.2378944 2.3194210
42 -18.5769239 12.2378944
43 3.6774620 -18.5769239
44 -13.4639077 3.6774620
45 -1.4073105 -13.4639077
46 13.9433930 -1.4073105
47 1.9175530 13.9433930
48 1.7214069 1.9175530
49 -14.7140856 1.7214069
50 6.7274408 -14.7140856
51 7.0979219 6.7274408
52 -15.3308004 7.0979219
53 4.4339402 -15.3308004
54 2.4048990 4.4339402
55 4.3458867 2.4048990
56 5.0417024 4.3458867
57 -11.9387601 5.0417024
> 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/7958q1258720956.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/8y46c1258720956.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/991291258720956.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/101d0o1258720956.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/11f5z11258720956.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/120cjn1258720956.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/13dwx71258720956.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/1441wn1258720956.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/15tguk1258720956.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/16h8fs1258720956.tab")
+ }
>
> system("convert tmp/1sa701258720956.ps tmp/1sa701258720956.png")
> system("convert tmp/2f8s81258720956.ps tmp/2f8s81258720956.png")
> system("convert tmp/36scx1258720956.ps tmp/36scx1258720956.png")
> system("convert tmp/47qdz1258720956.ps tmp/47qdz1258720956.png")
> system("convert tmp/5i64d1258720956.ps tmp/5i64d1258720956.png")
> system("convert tmp/6g06e1258720956.ps tmp/6g06e1258720956.png")
> system("convert tmp/7958q1258720956.ps tmp/7958q1258720956.png")
> system("convert tmp/8y46c1258720956.ps tmp/8y46c1258720956.png")
> system("convert tmp/991291258720956.ps tmp/991291258720956.png")
> system("convert tmp/101d0o1258720956.ps tmp/101d0o1258720956.png")
>
>
> proc.time()
user system elapsed
2.389 1.567 3.165