R version 2.12.0 (2010-10-15)
Copyright (C) 2010 The R Foundation for Statistical Computing
ISBN 3-900051-07-0
Platform: i486-pc-linux-gnu (32-bit)
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(153452
+ ,0
+ ,169422
+ ,174000
+ ,80900
+ ,35600
+ ,36700
+ ,173570
+ ,0
+ ,153452
+ ,169422
+ ,174000
+ ,80900
+ ,35600
+ ,193036
+ ,0
+ ,173570
+ ,153452
+ ,169422
+ ,174000
+ ,80900
+ ,174652
+ ,0
+ ,193036
+ ,173570
+ ,153452
+ ,169422
+ ,174000
+ ,105367
+ ,0
+ ,174652
+ ,193036
+ ,173570
+ ,153452
+ ,169422
+ ,95963
+ ,0
+ ,105367
+ ,174652
+ ,193036
+ ,173570
+ ,153452
+ ,82896
+ ,0
+ ,95963
+ ,105367
+ ,174652
+ ,193036
+ ,173570
+ ,121747
+ ,0
+ ,82896
+ ,95963
+ ,105367
+ ,174652
+ ,193036
+ ,120196
+ ,0
+ ,121747
+ ,82896
+ ,95963
+ ,105367
+ ,174652
+ ,103983
+ ,0
+ ,120196
+ ,121747
+ ,82896
+ ,95963
+ ,105367
+ ,81103
+ ,0
+ ,103983
+ ,120196
+ ,121747
+ ,82896
+ ,95963
+ ,70944
+ ,0
+ ,81103
+ ,103983
+ ,120196
+ ,121747
+ ,82896
+ ,57248
+ ,0
+ ,70944
+ ,81103
+ ,103983
+ ,120196
+ ,121747
+ ,47830
+ ,0
+ ,57248
+ ,70944
+ ,81103
+ ,103983
+ ,120196
+ ,60095
+ ,0
+ ,47830
+ ,57248
+ ,70944
+ ,81103
+ ,103983
+ ,60931
+ ,0
+ ,60095
+ ,47830
+ ,57248
+ ,70944
+ ,81103
+ ,82955
+ ,0
+ ,60931
+ ,60095
+ ,47830
+ ,57248
+ ,70944
+ ,99559
+ ,0
+ ,82955
+ ,60931
+ ,60095
+ ,47830
+ ,57248
+ ,77911
+ ,0
+ ,99559
+ ,82955
+ ,60931
+ ,60095
+ ,47830
+ ,70753
+ ,0
+ ,77911
+ ,99559
+ ,82955
+ ,60931
+ ,60095
+ ,69287
+ ,0
+ ,70753
+ ,77911
+ ,99559
+ ,82955
+ ,60931
+ ,88426
+ ,0
+ ,69287
+ ,70753
+ ,77911
+ ,99559
+ ,82955
+ ,91756
+ ,1
+ ,88426
+ ,69287
+ ,70753
+ ,77911
+ ,99559
+ ,96933
+ ,1
+ ,91756
+ ,88426
+ ,69287
+ ,70753
+ ,77911
+ ,174484
+ ,1
+ ,96933
+ ,91756
+ ,88426
+ ,69287
+ ,70753
+ ,232595
+ ,1
+ ,174484
+ ,96933
+ ,91756
+ ,88426
+ ,69287
+ ,266197
+ ,1
+ ,232595
+ ,174484
+ ,96933
+ ,91756
+ ,88426
+ ,290435
+ ,1
+ ,266197
+ ,232595
+ ,174484
+ ,96933
+ ,91756
+ ,304296
+ ,1
+ ,290435
+ ,266197
+ ,232595
+ ,174484
+ ,96933
+ ,322310
+ ,1
+ ,304296
+ ,290435
+ ,266197
+ ,232595
+ ,174484
+ ,415555
+ ,1
+ ,322310
+ ,304296
+ ,290435
+ ,266197
+ ,232595
+ ,490042
+ ,1
+ ,415555
+ ,322310
+ ,304296
+ ,290435
+ ,266197
+ ,545109
+ ,1
+ ,490042
+ ,415555
+ ,322310
+ ,304296
+ ,290435
+ ,545720
+ ,1
+ ,545109
+ ,490042
+ ,415555
+ ,322310
+ ,304296
+ ,505944
+ ,1
+ ,545720
+ ,545109
+ ,490042
+ ,415555
+ ,322310
+ ,477930
+ ,1
+ ,505944
+ ,545720
+ ,545109
+ ,490042
+ ,415555
+ ,466106
+ ,1
+ ,477930
+ ,505944
+ ,545720
+ ,545109
+ ,490042
+ ,424476
+ ,1
+ ,466106
+ ,477930
+ ,505944
+ ,545720
+ ,545109
+ ,383018
+ ,1
+ ,424476
+ ,466106
+ ,477930
+ ,505944
+ ,545720
+ ,364696
+ ,1
+ ,383018
+ ,424476
+ ,466106
+ ,477930
+ ,505944
+ ,391116
+ ,1
+ ,364696
+ ,383018
+ ,424476
+ ,466106
+ ,477930
+ ,435721
+ ,1
+ ,391116
+ ,364696
+ ,383018
+ ,424476
+ ,466106
+ ,511435
+ ,1
+ ,435721
+ ,391116
+ ,364696
+ ,383018
+ ,424476
+ ,553997
+ ,1
+ ,511435
+ ,435721
+ ,391116
+ ,364696
+ ,383018
+ ,555252
+ ,1
+ ,553997
+ ,511435
+ ,435721
+ ,391116
+ ,364696
+ ,544897
+ ,1
+ ,555252
+ ,553997
+ ,511435
+ ,435721
+ ,391116
+ ,540562
+ ,1
+ ,544897
+ ,555252
+ ,553997
+ ,511435
+ ,435721
+ ,505282
+ ,1
+ ,540562
+ ,544897
+ ,555252
+ ,553997
+ ,511435
+ ,507626
+ ,1
+ ,505282
+ ,540562
+ ,544897
+ ,555252
+ ,553997
+ ,474427
+ ,1
+ ,507626
+ ,505282
+ ,540562
+ ,544897
+ ,555252
+ ,469740
+ ,1
+ ,474427
+ ,507626
+ ,505282
+ ,540562
+ ,544897
+ ,491480
+ ,1
+ ,469740
+ ,474427
+ ,507626
+ ,505282
+ ,540562
+ ,538974
+ ,1
+ ,491480
+ ,469740
+ ,474427
+ ,507626
+ ,505282
+ ,576612
+ ,1
+ ,538974
+ ,491480
+ ,469740
+ ,474427
+ ,507626)
+ ,dim=c(7
+ ,54)
+ ,dimnames=list(c('Werklozen'
+ ,'Oliecrisis'
+ ,'Y1'
+ ,'Y2'
+ ,'Y3'
+ ,'Y4'
+ ,'Y5')
+ ,1:54))
> y <- array(NA,dim=c(7,54),dimnames=list(c('Werklozen','Oliecrisis','Y1','Y2','Y3','Y4','Y5'),1:54))
> 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 = 'Do not include Seasonal 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
> 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
Werklozen Oliecrisis Y1 Y2 Y3 Y4 Y5 t
1 153452 0 169422 174000 80900 35600 36700 1
2 173570 0 153452 169422 174000 80900 35600 2
3 193036 0 173570 153452 169422 174000 80900 3
4 174652 0 193036 173570 153452 169422 174000 4
5 105367 0 174652 193036 173570 153452 169422 5
6 95963 0 105367 174652 193036 173570 153452 6
7 82896 0 95963 105367 174652 193036 173570 7
8 121747 0 82896 95963 105367 174652 193036 8
9 120196 0 121747 82896 95963 105367 174652 9
10 103983 0 120196 121747 82896 95963 105367 10
11 81103 0 103983 120196 121747 82896 95963 11
12 70944 0 81103 103983 120196 121747 82896 12
13 57248 0 70944 81103 103983 120196 121747 13
14 47830 0 57248 70944 81103 103983 120196 14
15 60095 0 47830 57248 70944 81103 103983 15
16 60931 0 60095 47830 57248 70944 81103 16
17 82955 0 60931 60095 47830 57248 70944 17
18 99559 0 82955 60931 60095 47830 57248 18
19 77911 0 99559 82955 60931 60095 47830 19
20 70753 0 77911 99559 82955 60931 60095 20
21 69287 0 70753 77911 99559 82955 60931 21
22 88426 0 69287 70753 77911 99559 82955 22
23 91756 1 88426 69287 70753 77911 99559 23
24 96933 1 91756 88426 69287 70753 77911 24
25 174484 1 96933 91756 88426 69287 70753 25
26 232595 1 174484 96933 91756 88426 69287 26
27 266197 1 232595 174484 96933 91756 88426 27
28 290435 1 266197 232595 174484 96933 91756 28
29 304296 1 290435 266197 232595 174484 96933 29
30 322310 1 304296 290435 266197 232595 174484 30
31 415555 1 322310 304296 290435 266197 232595 31
32 490042 1 415555 322310 304296 290435 266197 32
33 545109 1 490042 415555 322310 304296 290435 33
34 545720 1 545109 490042 415555 322310 304296 34
35 505944 1 545720 545109 490042 415555 322310 35
36 477930 1 505944 545720 545109 490042 415555 36
37 466106 1 477930 505944 545720 545109 490042 37
38 424476 1 466106 477930 505944 545720 545109 38
39 383018 1 424476 466106 477930 505944 545720 39
40 364696 1 383018 424476 466106 477930 505944 40
41 391116 1 364696 383018 424476 466106 477930 41
42 435721 1 391116 364696 383018 424476 466106 42
43 511435 1 435721 391116 364696 383018 424476 43
44 553997 1 511435 435721 391116 364696 383018 44
45 555252 1 553997 511435 435721 391116 364696 45
46 544897 1 555252 553997 511435 435721 391116 46
47 540562 1 544897 555252 553997 511435 435721 47
48 505282 1 540562 544897 555252 553997 511435 48
49 507626 1 505282 540562 544897 555252 553997 49
50 474427 1 507626 505282 540562 544897 555252 50
51 469740 1 474427 507626 505282 540562 544897 51
52 491480 1 469740 474427 507626 505282 540562 52
53 538974 1 491480 469740 474427 507626 505282 53
54 576612 1 538974 491480 469740 474427 507626 54
> k <- length(x[1,])
> df <- as.data.frame(x)
> (mylm <- lm(df))
Call:
lm(formula = df)
Coefficients:
(Intercept) Oliecrisis Y1 Y2 Y3 Y4
2.989e+03 2.391e+04 1.499e+00 -6.735e-01 -2.213e-02 1.671e-01
Y5 t
-8.113e-02 6.595e+02
> (mysum <- summary(mylm))
Call:
lm(formula = df)
Residuals:
Min 1Q Median 3Q Max
-40739 -14135 -1684 14488 70878
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 2.989e+03 7.674e+03 0.389 0.6988
Oliecrisis 2.391e+04 1.525e+04 1.568 0.1238
Y1 1.499e+00 1.436e-01 10.437 1.03e-13 ***
Y2 -6.735e-01 2.446e-01 -2.754 0.0084 **
Y3 -2.213e-02 2.548e-01 -0.087 0.9312
Y4 1.671e-01 2.356e-01 0.710 0.4815
Y5 -8.113e-02 1.282e-01 -0.633 0.5300
t 6.595e+02 6.172e+02 1.069 0.2908
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 25900 on 46 degrees of freedom
Multiple R-squared: 0.9842, Adjusted R-squared: 0.9818
F-statistic: 409.1 on 7 and 46 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.74335312 0.5132938 0.2566469
[2,] 0.61782586 0.7643483 0.3821741
[3,] 0.48744804 0.9748961 0.5125520
[4,] 0.35341222 0.7068244 0.6465878
[5,] 0.27901558 0.5580312 0.7209844
[6,] 0.18494771 0.3698954 0.8150523
[7,] 0.23437275 0.4687455 0.7656273
[8,] 0.22742667 0.4548533 0.7725733
[9,] 0.17729703 0.3545941 0.8227030
[10,] 0.15934202 0.3186840 0.8406580
[11,] 0.11236896 0.2247379 0.8876310
[12,] 0.10447493 0.2089499 0.8955251
[13,] 0.09597969 0.1919594 0.9040203
[14,] 0.09072071 0.1814414 0.9092793
[15,] 0.36219578 0.7243916 0.6378042
[16,] 0.36842653 0.7368531 0.6315735
[17,] 0.41804399 0.8360880 0.5819560
[18,] 0.39463605 0.7892721 0.6053639
[19,] 0.49627289 0.9925458 0.5037271
[20,] 0.84044224 0.3191155 0.1595578
[21,] 0.92823100 0.1435380 0.0717690
[22,] 0.89711690 0.2057662 0.1028831
[23,] 0.85920475 0.2815905 0.1407953
[24,] 0.83326836 0.3334633 0.1667316
[25,] 0.86376591 0.2724682 0.1362341
[26,] 0.83230052 0.3353990 0.1676995
[27,] 0.90698540 0.1860292 0.0930146
[28,] 0.87810283 0.2437943 0.1218972
[29,] 0.80040543 0.3991891 0.1995946
[30,] 0.72361046 0.5527791 0.2763895
[31,] 0.63464266 0.7307147 0.3653573
[32,] 0.54668525 0.9066295 0.4533147
[33,] 0.47915419 0.9583084 0.5208458
> postscript(file="/var/www/rcomp/tmp/10pi41292687632.ps",horizontal=F,onefile=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/rcomp/tmp/20pi41292687632.ps",horizontal=F,onefile=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/rcomp/tmp/3tyz71292687632.ps",horizontal=F,onefile=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/rcomp/tmp/4tyz71292687632.ps",horizontal=F,onefile=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/rcomp/tmp/5tyz71292687632.ps",horizontal=F,onefile=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 = 54
Frequency = 1
1 2 3 4 5 6
11880.99264 46591.54868 12501.68170 -14202.66137 -40739.34959 36432.85999
7 8 9 10 11 12
-11891.41169 42670.13652 -16690.49369 -9410.48849 -7413.50068 -2447.18861
13 14 15 16 17 18
-13933.68710 -8247.97443 10532.91124 -14478.11398 15150.52039 -617.39306
19 20 21 22 23 24
-35773.44465 1381.06323 -7842.11656 6546.10951 -39561.83749 -27737.48375
25 26 27 28 29 30
43725.25784 -14815.11072 -15627.76033 -2152.77886 -13904.43541 -3677.77387
31 32 33 34 35 36
70878.39031 16063.54249 21679.03485 -10559.44883 -27298.00323 390.51984
37 38 39 40 41 42
-42.20907 -39992.16610 -21599.79257 -5288.35432 18792.99413 15881.21550
43 44 45 46 47 48
45021.81273 3768.74114 -13349.32189 -1215.22485 2061.30383 -35298.41955
49 50 51 52 53 54
19358.72919 -40037.85509 5057.01194 16399.53953 23504.15006 9575.73254
> postscript(file="/var/www/rcomp/tmp/6mpya1292687632.ps",horizontal=F,onefile=F,pagecentre=F,paper="special",width=8.3333333333333,height=5.5555555555556)
> dum <- cbind(lag(myerror,k=1),myerror)
> dum
Time Series:
Start = 0
End = 54
Frequency = 1
lag(myerror, k = 1) myerror
0 11880.99264 NA
1 46591.54868 11880.99264
2 12501.68170 46591.54868
3 -14202.66137 12501.68170
4 -40739.34959 -14202.66137
5 36432.85999 -40739.34959
6 -11891.41169 36432.85999
7 42670.13652 -11891.41169
8 -16690.49369 42670.13652
9 -9410.48849 -16690.49369
10 -7413.50068 -9410.48849
11 -2447.18861 -7413.50068
12 -13933.68710 -2447.18861
13 -8247.97443 -13933.68710
14 10532.91124 -8247.97443
15 -14478.11398 10532.91124
16 15150.52039 -14478.11398
17 -617.39306 15150.52039
18 -35773.44465 -617.39306
19 1381.06323 -35773.44465
20 -7842.11656 1381.06323
21 6546.10951 -7842.11656
22 -39561.83749 6546.10951
23 -27737.48375 -39561.83749
24 43725.25784 -27737.48375
25 -14815.11072 43725.25784
26 -15627.76033 -14815.11072
27 -2152.77886 -15627.76033
28 -13904.43541 -2152.77886
29 -3677.77387 -13904.43541
30 70878.39031 -3677.77387
31 16063.54249 70878.39031
32 21679.03485 16063.54249
33 -10559.44883 21679.03485
34 -27298.00323 -10559.44883
35 390.51984 -27298.00323
36 -42.20907 390.51984
37 -39992.16610 -42.20907
38 -21599.79257 -39992.16610
39 -5288.35432 -21599.79257
40 18792.99413 -5288.35432
41 15881.21550 18792.99413
42 45021.81273 15881.21550
43 3768.74114 45021.81273
44 -13349.32189 3768.74114
45 -1215.22485 -13349.32189
46 2061.30383 -1215.22485
47 -35298.41955 2061.30383
48 19358.72919 -35298.41955
49 -40037.85509 19358.72919
50 5057.01194 -40037.85509
51 16399.53953 5057.01194
52 23504.15006 16399.53953
53 9575.73254 23504.15006
54 NA 9575.73254
> dum1 <- dum[2:length(myerror),]
> dum1
lag(myerror, k = 1) myerror
[1,] 46591.54868 11880.99264
[2,] 12501.68170 46591.54868
[3,] -14202.66137 12501.68170
[4,] -40739.34959 -14202.66137
[5,] 36432.85999 -40739.34959
[6,] -11891.41169 36432.85999
[7,] 42670.13652 -11891.41169
[8,] -16690.49369 42670.13652
[9,] -9410.48849 -16690.49369
[10,] -7413.50068 -9410.48849
[11,] -2447.18861 -7413.50068
[12,] -13933.68710 -2447.18861
[13,] -8247.97443 -13933.68710
[14,] 10532.91124 -8247.97443
[15,] -14478.11398 10532.91124
[16,] 15150.52039 -14478.11398
[17,] -617.39306 15150.52039
[18,] -35773.44465 -617.39306
[19,] 1381.06323 -35773.44465
[20,] -7842.11656 1381.06323
[21,] 6546.10951 -7842.11656
[22,] -39561.83749 6546.10951
[23,] -27737.48375 -39561.83749
[24,] 43725.25784 -27737.48375
[25,] -14815.11072 43725.25784
[26,] -15627.76033 -14815.11072
[27,] -2152.77886 -15627.76033
[28,] -13904.43541 -2152.77886
[29,] -3677.77387 -13904.43541
[30,] 70878.39031 -3677.77387
[31,] 16063.54249 70878.39031
[32,] 21679.03485 16063.54249
[33,] -10559.44883 21679.03485
[34,] -27298.00323 -10559.44883
[35,] 390.51984 -27298.00323
[36,] -42.20907 390.51984
[37,] -39992.16610 -42.20907
[38,] -21599.79257 -39992.16610
[39,] -5288.35432 -21599.79257
[40,] 18792.99413 -5288.35432
[41,] 15881.21550 18792.99413
[42,] 45021.81273 15881.21550
[43,] 3768.74114 45021.81273
[44,] -13349.32189 3768.74114
[45,] -1215.22485 -13349.32189
[46,] 2061.30383 -1215.22485
[47,] -35298.41955 2061.30383
[48,] 19358.72919 -35298.41955
[49,] -40037.85509 19358.72919
[50,] 5057.01194 -40037.85509
[51,] 16399.53953 5057.01194
[52,] 23504.15006 16399.53953
[53,] 9575.73254 23504.15006
> z <- as.data.frame(dum1)
> z
lag(myerror, k = 1) myerror
1 46591.54868 11880.99264
2 12501.68170 46591.54868
3 -14202.66137 12501.68170
4 -40739.34959 -14202.66137
5 36432.85999 -40739.34959
6 -11891.41169 36432.85999
7 42670.13652 -11891.41169
8 -16690.49369 42670.13652
9 -9410.48849 -16690.49369
10 -7413.50068 -9410.48849
11 -2447.18861 -7413.50068
12 -13933.68710 -2447.18861
13 -8247.97443 -13933.68710
14 10532.91124 -8247.97443
15 -14478.11398 10532.91124
16 15150.52039 -14478.11398
17 -617.39306 15150.52039
18 -35773.44465 -617.39306
19 1381.06323 -35773.44465
20 -7842.11656 1381.06323
21 6546.10951 -7842.11656
22 -39561.83749 6546.10951
23 -27737.48375 -39561.83749
24 43725.25784 -27737.48375
25 -14815.11072 43725.25784
26 -15627.76033 -14815.11072
27 -2152.77886 -15627.76033
28 -13904.43541 -2152.77886
29 -3677.77387 -13904.43541
30 70878.39031 -3677.77387
31 16063.54249 70878.39031
32 21679.03485 16063.54249
33 -10559.44883 21679.03485
34 -27298.00323 -10559.44883
35 390.51984 -27298.00323
36 -42.20907 390.51984
37 -39992.16610 -42.20907
38 -21599.79257 -39992.16610
39 -5288.35432 -21599.79257
40 18792.99413 -5288.35432
41 15881.21550 18792.99413
42 45021.81273 15881.21550
43 3768.74114 45021.81273
44 -13349.32189 3768.74114
45 -1215.22485 -13349.32189
46 2061.30383 -1215.22485
47 -35298.41955 2061.30383
48 19358.72919 -35298.41955
49 -40037.85509 19358.72919
50 5057.01194 -40037.85509
51 16399.53953 5057.01194
52 23504.15006 16399.53953
53 9575.73254 23504.15006
> 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/rcomp/tmp/7wgxd1292687632.ps",horizontal=F,onefile=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/rcomp/tmp/8wgxd1292687632.ps",horizontal=F,onefile=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/rcomp/tmp/9wgxd1292687632.ps",horizontal=F,onefile=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/rcomp/tmp/107pxy1292687632.ps",horizontal=F,onefile=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/rcomp/createtable file can be downloaded at http://www.wessa.net/cretab
> load(file="/var/www/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/rcomp/tmp/11yt2s1292687632.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/rcomp/tmp/12wrca1292687632.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/rcomp/tmp/13a0ri1292687632.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/rcomp/tmp/14vj8o1292687632.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/rcomp/tmp/15hjpu1292687632.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/rcomp/tmp/16k25i1292687632.tab")
+ }
> try(system("convert tmp/10pi41292687632.ps tmp/10pi41292687632.png",intern=TRUE))
character(0)
> try(system("convert tmp/20pi41292687632.ps tmp/20pi41292687632.png",intern=TRUE))
character(0)
> try(system("convert tmp/3tyz71292687632.ps tmp/3tyz71292687632.png",intern=TRUE))
character(0)
> try(system("convert tmp/4tyz71292687632.ps tmp/4tyz71292687632.png",intern=TRUE))
character(0)
> try(system("convert tmp/5tyz71292687632.ps tmp/5tyz71292687632.png",intern=TRUE))
character(0)
> try(system("convert tmp/6mpya1292687632.ps tmp/6mpya1292687632.png",intern=TRUE))
character(0)
> try(system("convert tmp/7wgxd1292687632.ps tmp/7wgxd1292687632.png",intern=TRUE))
character(0)
> try(system("convert tmp/8wgxd1292687632.ps tmp/8wgxd1292687632.png",intern=TRUE))
character(0)
> try(system("convert tmp/9wgxd1292687632.ps tmp/9wgxd1292687632.png",intern=TRUE))
character(0)
> try(system("convert tmp/107pxy1292687632.ps tmp/107pxy1292687632.png",intern=TRUE))
character(0)
>
>
> proc.time()
user system elapsed
3.01 0.89 3.87