R version 2.13.0 (2011-04-13)
Copyright (C) 2011 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(118.49
+ ,548
+ ,4.2
+ ,23118
+ ,2075
+ ,118.31
+ ,563
+ ,3.88
+ ,22849
+ ,2294
+ ,117.99
+ ,581
+ ,4.11
+ ,20198
+ ,2670
+ ,118.09
+ ,572
+ ,4.22
+ ,18130
+ ,2242
+ ,117.95
+ ,519
+ ,4.14
+ ,25597
+ ,2764
+ ,117.59
+ ,521
+ ,4.21
+ ,28785
+ ,2409
+ ,117.2
+ ,531
+ ,4.29
+ ,28229
+ ,2187
+ ,116.91
+ ,540
+ ,4.21
+ ,33474
+ ,2221
+ ,116.33
+ ,548
+ ,4.21
+ ,28287
+ ,1777
+ ,115.66
+ ,556
+ ,4.14
+ ,27297
+ ,1956
+ ,115
+ ,551
+ ,3.99
+ ,16167
+ ,1493
+ ,114.55
+ ,549
+ ,3.48
+ ,22380
+ ,1719
+ ,114.41
+ ,564
+ ,3.21
+ ,24256
+ ,3471
+ ,114.25
+ ,586
+ ,3.12
+ ,19573
+ ,4894
+ ,113.89
+ ,604
+ ,3.03
+ ,21553
+ ,4242
+ ,113.82
+ ,601
+ ,3.29
+ ,21359
+ ,3595
+ ,113.77
+ ,545
+ ,3.47
+ ,29586
+ ,3762
+ ,113.78
+ ,537
+ ,3.31
+ ,27186
+ ,3055
+ ,113.33
+ ,552
+ ,3.54
+ ,32066
+ ,2503
+ ,112.94
+ ,563
+ ,3.63
+ ,34670
+ ,2327
+ ,112.52
+ ,575
+ ,3.73
+ ,25985
+ ,2389
+ ,112.05
+ ,580
+ ,3.75
+ ,27561
+ ,2923
+ ,111.54
+ ,575
+ ,3.61
+ ,14538
+ ,2624
+ ,111.36
+ ,558
+ ,3.64
+ ,18730
+ ,2424
+ ,111.07
+ ,564
+ ,3.68
+ ,22485
+ ,2592
+ ,111.02
+ ,581
+ ,3.72
+ ,20036
+ ,2859
+ ,111.31
+ ,597
+ ,3.77
+ ,16971
+ ,2349
+ ,110.97
+ ,587
+ ,3.92
+ ,19028
+ ,2524
+ ,111.04
+ ,536
+ ,4.12
+ ,22759
+ ,2622
+ ,111.25
+ ,524
+ ,4.03
+ ,20516
+ ,2362
+ ,111.33
+ ,537
+ ,3.93
+ ,26195
+ ,2251
+ ,111.1
+ ,536
+ ,4.03
+ ,27786
+ ,3071
+ ,111.74
+ ,533
+ ,4.24
+ ,24090
+ ,2859
+ ,111.36
+ ,528
+ ,4.13
+ ,25447
+ ,2645
+ ,111.25
+ ,516
+ ,3.87
+ ,11509
+ ,3133
+ ,111.49
+ ,502
+ ,4.26
+ ,15572
+ ,2575
+ ,112.16
+ ,506
+ ,4.46
+ ,22518
+ ,2583
+ ,112.36
+ ,518
+ ,4.56
+ ,20520
+ ,3200
+ ,112.18
+ ,534
+ ,4.58
+ ,17789
+ ,2875
+ ,112.87
+ ,528
+ ,4.85
+ ,20205
+ ,3014
+ ,112.28
+ ,478
+ ,4.84
+ ,26835
+ ,2925
+ ,111.66
+ ,469
+ ,4.51
+ ,25826
+ ,3373
+ ,110.67
+ ,490
+ ,4.37
+ ,31934
+ ,2925
+ ,110.42
+ ,493
+ ,4.23
+ ,30019
+ ,2591
+ ,109.62
+ ,508
+ ,4.23
+ ,30111
+ ,2814
+ ,108.84
+ ,517
+ ,4.25
+ ,31566
+ ,3641
+ ,108.4
+ ,514
+ ,4.41
+ ,12738
+ ,2578
+ ,108.1
+ ,510
+ ,4.28
+ ,19814
+ ,3129
+ ,107.1
+ ,527
+ ,4.42
+ ,24776
+ ,2849
+ ,106.54
+ ,542
+ ,4.39
+ ,20424
+ ,3534
+ ,106.44
+ ,565
+ ,4.44
+ ,18688
+ ,2617
+ ,106.57
+ ,555
+ ,4.62
+ ,20418
+ ,3016
+ ,106.12
+ ,499
+ ,4.64
+ ,25778
+ ,3483
+ ,106.13
+ ,511
+ ,4.34
+ ,25100
+ ,3014
+ ,106.26
+ ,526
+ ,4.22
+ ,25859
+ ,3179
+ ,105.78
+ ,532
+ ,4.01
+ ,30651
+ ,2907
+ ,105.77
+ ,549
+ ,4.11
+ ,26551
+ ,2770
+ ,105.2
+ ,561
+ ,4.06
+ ,31124
+ ,3498
+ ,105.15
+ ,557
+ ,3.82
+ ,9367
+ ,3417
+ ,105.01
+ ,566
+ ,3.76
+ ,17382
+ ,3324)
+ ,dim=c(5
+ ,60)
+ ,dimnames=list(c('CPI'
+ ,'werkloosheid'
+ ,'OLO'
+ ,'voertuigen'
+ ,'bouw')
+ ,1:60))
> y <- array(NA,dim=c(5,60),dimnames=list(c('CPI','werkloosheid','OLO','voertuigen','bouw'),1:60))
> 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
CPI werkloosheid OLO voertuigen bouw t
1 118.49 548 4.20 23118 2075 1
2 118.31 563 3.88 22849 2294 2
3 117.99 581 4.11 20198 2670 3
4 118.09 572 4.22 18130 2242 4
5 117.95 519 4.14 25597 2764 5
6 117.59 521 4.21 28785 2409 6
7 117.20 531 4.29 28229 2187 7
8 116.91 540 4.21 33474 2221 8
9 116.33 548 4.21 28287 1777 9
10 115.66 556 4.14 27297 1956 10
11 115.00 551 3.99 16167 1493 11
12 114.55 549 3.48 22380 1719 12
13 114.41 564 3.21 24256 3471 13
14 114.25 586 3.12 19573 4894 14
15 113.89 604 3.03 21553 4242 15
16 113.82 601 3.29 21359 3595 16
17 113.77 545 3.47 29586 3762 17
18 113.78 537 3.31 27186 3055 18
19 113.33 552 3.54 32066 2503 19
20 112.94 563 3.63 34670 2327 20
21 112.52 575 3.73 25985 2389 21
22 112.05 580 3.75 27561 2923 22
23 111.54 575 3.61 14538 2624 23
24 111.36 558 3.64 18730 2424 24
25 111.07 564 3.68 22485 2592 25
26 111.02 581 3.72 20036 2859 26
27 111.31 597 3.77 16971 2349 27
28 110.97 587 3.92 19028 2524 28
29 111.04 536 4.12 22759 2622 29
30 111.25 524 4.03 20516 2362 30
31 111.33 537 3.93 26195 2251 31
32 111.10 536 4.03 27786 3071 32
33 111.74 533 4.24 24090 2859 33
34 111.36 528 4.13 25447 2645 34
35 111.25 516 3.87 11509 3133 35
36 111.49 502 4.26 15572 2575 36
37 112.16 506 4.46 22518 2583 37
38 112.36 518 4.56 20520 3200 38
39 112.18 534 4.58 17789 2875 39
40 112.87 528 4.85 20205 3014 40
41 112.28 478 4.84 26835 2925 41
42 111.66 469 4.51 25826 3373 42
43 110.67 490 4.37 31934 2925 43
44 110.42 493 4.23 30019 2591 44
45 109.62 508 4.23 30111 2814 45
46 108.84 517 4.25 31566 3641 46
47 108.40 514 4.41 12738 2578 47
48 108.10 510 4.28 19814 3129 48
49 107.10 527 4.42 24776 2849 49
50 106.54 542 4.39 20424 3534 50
51 106.44 565 4.44 18688 2617 51
52 106.57 555 4.62 20418 3016 52
53 106.12 499 4.64 25778 3483 53
54 106.13 511 4.34 25100 3014 54
55 106.26 526 4.22 25859 3179 55
56 105.78 532 4.01 30651 2907 56
57 105.77 549 4.11 26551 2770 57
58 105.20 561 4.06 31124 3498 58
59 105.15 557 3.82 9367 3417 59
60 105.01 566 3.76 17382 3324 60
> k <- length(x[1,])
> df <- as.data.frame(x)
> (mylm <- lm(df))
Call:
lm(formula = df)
Coefficients:
(Intercept) werkloosheid OLO voertuigen bouw
1.207e+02 -1.832e-02 1.737e+00 -1.468e-05 5.562e-04
t
-2.424e-01
> (mysum <- summary(mylm))
Call:
lm(formula = df)
Residuals:
Min 1Q Median 3Q Max
-2.22471 -0.44604 0.01431 0.48009 1.71929
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 1.207e+02 3.695e+00 32.666 < 2e-16 ***
werkloosheid -1.832e-02 4.469e-03 -4.099 0.000141 ***
OLO 1.737e+00 3.473e-01 5.000 6.41e-06 ***
voertuigen -1.468e-05 1.786e-05 -0.822 0.414793
bouw 5.562e-04 1.897e-04 2.932 0.004933 **
t -2.424e-01 7.129e-03 -34.004 < 2e-16 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 0.7209 on 54 degrees of freedom
Multiple R-squared: 0.9663, Adjusted R-squared: 0.9632
F-statistic: 309.5 on 5 and 54 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,] 2.037781e-02 4.075562e-02 0.97962219
[2,] 2.326005e-02 4.652010e-02 0.97673995
[3,] 7.676110e-03 1.535222e-02 0.99232389
[4,] 2.101218e-03 4.202437e-03 0.99789878
[5,] 7.043497e-04 1.408699e-03 0.99929565
[6,] 1.826491e-04 3.652983e-04 0.99981735
[7,] 7.448953e-05 1.489791e-04 0.99992551
[8,] 5.397457e-05 1.079491e-04 0.99994603
[9,] 1.312392e-05 2.624784e-05 0.99998688
[10,] 9.124199e-05 1.824840e-04 0.99990876
[11,] 4.186180e-05 8.372360e-05 0.99995814
[12,] 1.285836e-05 2.571672e-05 0.99998714
[13,] 3.722661e-06 7.445322e-06 0.99999628
[14,] 1.639218e-06 3.278437e-06 0.99999836
[15,] 4.712718e-07 9.425435e-07 0.99999953
[16,] 1.421402e-07 2.842803e-07 0.99999986
[17,] 4.851242e-08 9.702485e-08 0.99999995
[18,] 1.681121e-08 3.362243e-08 0.99999998
[19,] 4.060009e-07 8.120017e-07 0.99999959
[20,] 2.232259e-07 4.464519e-07 0.99999978
[21,] 2.071881e-07 4.143762e-07 0.99999979
[22,] 1.245282e-06 2.490565e-06 0.99999875
[23,] 1.611577e-05 3.223153e-05 0.99998388
[24,] 2.226024e-05 4.452047e-05 0.99997774
[25,] 1.109818e-04 2.219636e-04 0.99988902
[26,] 4.231576e-04 8.463151e-04 0.99957684
[27,] 1.649664e-03 3.299328e-03 0.99835034
[28,] 4.788062e-03 9.576124e-03 0.99521194
[29,] 1.150540e-02 2.301080e-02 0.98849460
[30,] 1.238237e-02 2.476475e-02 0.98761763
[31,] 1.666662e-02 3.333324e-02 0.98333338
[32,] 1.749299e-01 3.498598e-01 0.82507012
[33,] 5.930448e-01 8.139103e-01 0.40695515
[34,] 8.739736e-01 2.520529e-01 0.12602643
[35,] 9.007214e-01 1.985573e-01 0.09927864
[36,] 9.287464e-01 1.425072e-01 0.07125361
[37,] 9.295442e-01 1.409117e-01 0.07045583
[38,] 9.610932e-01 7.781367e-02 0.03890684
[39,] 9.522945e-01 9.541109e-02 0.04770555
[40,] 9.936003e-01 1.279944e-02 0.00639972
[41,] 9.880824e-01 2.383517e-02 0.01191759
[42,] 9.749404e-01 5.011914e-02 0.02505957
[43,] 9.633306e-01 7.333888e-02 0.03666944
> postscript(file="/var/wessaorg/rcomp/tmp/100gz1321793415.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/wessaorg/rcomp/tmp/28vu51321793415.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/wessaorg/rcomp/tmp/38wpe1321793415.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/wessaorg/rcomp/tmp/468dw1321793415.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/wessaorg/rcomp/tmp/5mcic1321793415.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 = 60
Frequency = 1
1 2 3 4 5 6
-0.054997118 0.712196484 0.316824345 0.511050016 -0.399111849 -0.357388838
7 8 9 10 11 12
-0.345424493 -0.031134542 -0.051355747 -0.324910856 -0.479391280 0.127650260
13 14 15 16 17 18
0.026841143 -0.191651620 0.568485742 0.591405838 -0.526667276 0.215112960
19 20 21 22 23 24
0.261472615 0.295180842 0.001773072 -0.442830581 -0.583678644 -0.711977629
25 26 27 28 29 30
-0.757451505 -0.507561030 0.469766648 -0.138637773 -1.107473800 -0.606850384
31 32 33 34 35 36
0.272460579 -0.339838229 0.186577167 0.287406827 0.175600964 0.094233455
37 38 39 40 41 42
0.830066108 0.946116770 1.407554267 1.719294771 0.620058203 0.386786826
43 44 45 46 47 48
0.605837665 1.054024920 0.648520518 -0.197556815 -0.413075775 -0.520743476
49 50 51 52 53 54
-0.981520767 -1.417108749 -0.455686391 -0.775581380 -2.224706428 -0.980540211
55 56 57 58 59 60
-0.205581279 0.253087228 0.639241006 0.280515820 0.542222356 1.083069047
> postscript(file="/var/wessaorg/rcomp/tmp/6c7na1321793415.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 = 60
Frequency = 1
lag(myerror, k = 1) myerror
0 -0.054997118 NA
1 0.712196484 -0.054997118
2 0.316824345 0.712196484
3 0.511050016 0.316824345
4 -0.399111849 0.511050016
5 -0.357388838 -0.399111849
6 -0.345424493 -0.357388838
7 -0.031134542 -0.345424493
8 -0.051355747 -0.031134542
9 -0.324910856 -0.051355747
10 -0.479391280 -0.324910856
11 0.127650260 -0.479391280
12 0.026841143 0.127650260
13 -0.191651620 0.026841143
14 0.568485742 -0.191651620
15 0.591405838 0.568485742
16 -0.526667276 0.591405838
17 0.215112960 -0.526667276
18 0.261472615 0.215112960
19 0.295180842 0.261472615
20 0.001773072 0.295180842
21 -0.442830581 0.001773072
22 -0.583678644 -0.442830581
23 -0.711977629 -0.583678644
24 -0.757451505 -0.711977629
25 -0.507561030 -0.757451505
26 0.469766648 -0.507561030
27 -0.138637773 0.469766648
28 -1.107473800 -0.138637773
29 -0.606850384 -1.107473800
30 0.272460579 -0.606850384
31 -0.339838229 0.272460579
32 0.186577167 -0.339838229
33 0.287406827 0.186577167
34 0.175600964 0.287406827
35 0.094233455 0.175600964
36 0.830066108 0.094233455
37 0.946116770 0.830066108
38 1.407554267 0.946116770
39 1.719294771 1.407554267
40 0.620058203 1.719294771
41 0.386786826 0.620058203
42 0.605837665 0.386786826
43 1.054024920 0.605837665
44 0.648520518 1.054024920
45 -0.197556815 0.648520518
46 -0.413075775 -0.197556815
47 -0.520743476 -0.413075775
48 -0.981520767 -0.520743476
49 -1.417108749 -0.981520767
50 -0.455686391 -1.417108749
51 -0.775581380 -0.455686391
52 -2.224706428 -0.775581380
53 -0.980540211 -2.224706428
54 -0.205581279 -0.980540211
55 0.253087228 -0.205581279
56 0.639241006 0.253087228
57 0.280515820 0.639241006
58 0.542222356 0.280515820
59 1.083069047 0.542222356
60 NA 1.083069047
> dum1 <- dum[2:length(myerror),]
> dum1
lag(myerror, k = 1) myerror
[1,] 0.712196484 -0.054997118
[2,] 0.316824345 0.712196484
[3,] 0.511050016 0.316824345
[4,] -0.399111849 0.511050016
[5,] -0.357388838 -0.399111849
[6,] -0.345424493 -0.357388838
[7,] -0.031134542 -0.345424493
[8,] -0.051355747 -0.031134542
[9,] -0.324910856 -0.051355747
[10,] -0.479391280 -0.324910856
[11,] 0.127650260 -0.479391280
[12,] 0.026841143 0.127650260
[13,] -0.191651620 0.026841143
[14,] 0.568485742 -0.191651620
[15,] 0.591405838 0.568485742
[16,] -0.526667276 0.591405838
[17,] 0.215112960 -0.526667276
[18,] 0.261472615 0.215112960
[19,] 0.295180842 0.261472615
[20,] 0.001773072 0.295180842
[21,] -0.442830581 0.001773072
[22,] -0.583678644 -0.442830581
[23,] -0.711977629 -0.583678644
[24,] -0.757451505 -0.711977629
[25,] -0.507561030 -0.757451505
[26,] 0.469766648 -0.507561030
[27,] -0.138637773 0.469766648
[28,] -1.107473800 -0.138637773
[29,] -0.606850384 -1.107473800
[30,] 0.272460579 -0.606850384
[31,] -0.339838229 0.272460579
[32,] 0.186577167 -0.339838229
[33,] 0.287406827 0.186577167
[34,] 0.175600964 0.287406827
[35,] 0.094233455 0.175600964
[36,] 0.830066108 0.094233455
[37,] 0.946116770 0.830066108
[38,] 1.407554267 0.946116770
[39,] 1.719294771 1.407554267
[40,] 0.620058203 1.719294771
[41,] 0.386786826 0.620058203
[42,] 0.605837665 0.386786826
[43,] 1.054024920 0.605837665
[44,] 0.648520518 1.054024920
[45,] -0.197556815 0.648520518
[46,] -0.413075775 -0.197556815
[47,] -0.520743476 -0.413075775
[48,] -0.981520767 -0.520743476
[49,] -1.417108749 -0.981520767
[50,] -0.455686391 -1.417108749
[51,] -0.775581380 -0.455686391
[52,] -2.224706428 -0.775581380
[53,] -0.980540211 -2.224706428
[54,] -0.205581279 -0.980540211
[55,] 0.253087228 -0.205581279
[56,] 0.639241006 0.253087228
[57,] 0.280515820 0.639241006
[58,] 0.542222356 0.280515820
[59,] 1.083069047 0.542222356
> z <- as.data.frame(dum1)
> z
lag(myerror, k = 1) myerror
1 0.712196484 -0.054997118
2 0.316824345 0.712196484
3 0.511050016 0.316824345
4 -0.399111849 0.511050016
5 -0.357388838 -0.399111849
6 -0.345424493 -0.357388838
7 -0.031134542 -0.345424493
8 -0.051355747 -0.031134542
9 -0.324910856 -0.051355747
10 -0.479391280 -0.324910856
11 0.127650260 -0.479391280
12 0.026841143 0.127650260
13 -0.191651620 0.026841143
14 0.568485742 -0.191651620
15 0.591405838 0.568485742
16 -0.526667276 0.591405838
17 0.215112960 -0.526667276
18 0.261472615 0.215112960
19 0.295180842 0.261472615
20 0.001773072 0.295180842
21 -0.442830581 0.001773072
22 -0.583678644 -0.442830581
23 -0.711977629 -0.583678644
24 -0.757451505 -0.711977629
25 -0.507561030 -0.757451505
26 0.469766648 -0.507561030
27 -0.138637773 0.469766648
28 -1.107473800 -0.138637773
29 -0.606850384 -1.107473800
30 0.272460579 -0.606850384
31 -0.339838229 0.272460579
32 0.186577167 -0.339838229
33 0.287406827 0.186577167
34 0.175600964 0.287406827
35 0.094233455 0.175600964
36 0.830066108 0.094233455
37 0.946116770 0.830066108
38 1.407554267 0.946116770
39 1.719294771 1.407554267
40 0.620058203 1.719294771
41 0.386786826 0.620058203
42 0.605837665 0.386786826
43 1.054024920 0.605837665
44 0.648520518 1.054024920
45 -0.197556815 0.648520518
46 -0.413075775 -0.197556815
47 -0.520743476 -0.413075775
48 -0.981520767 -0.520743476
49 -1.417108749 -0.981520767
50 -0.455686391 -1.417108749
51 -0.775581380 -0.455686391
52 -2.224706428 -0.775581380
53 -0.980540211 -2.224706428
54 -0.205581279 -0.980540211
55 0.253087228 -0.205581279
56 0.639241006 0.253087228
57 0.280515820 0.639241006
58 0.542222356 0.280515820
59 1.083069047 0.542222356
> 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/wessaorg/rcomp/tmp/7tz8a1321793415.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/wessaorg/rcomp/tmp/8ch1d1321793415.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/wessaorg/rcomp/tmp/9ihdq1321793415.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/wessaorg/rcomp/tmp/10c7vl1321793415.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/wessaorg/rcomp/createtable file can be downloaded at http://www.wessa.net/cretab
> load(file="/var/wessaorg/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/wessaorg/rcomp/tmp/114l0a1321793415.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/wessaorg/rcomp/tmp/12ud7g1321793415.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/wessaorg/rcomp/tmp/13ppt21321793415.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/wessaorg/rcomp/tmp/14xgcq1321793415.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/wessaorg/rcomp/tmp/15fl1s1321793415.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/wessaorg/rcomp/tmp/16226l1321793415.tab")
+ }
>
> try(system("convert tmp/100gz1321793415.ps tmp/100gz1321793415.png",intern=TRUE))
character(0)
> try(system("convert tmp/28vu51321793415.ps tmp/28vu51321793415.png",intern=TRUE))
character(0)
> try(system("convert tmp/38wpe1321793415.ps tmp/38wpe1321793415.png",intern=TRUE))
character(0)
> try(system("convert tmp/468dw1321793415.ps tmp/468dw1321793415.png",intern=TRUE))
character(0)
> try(system("convert tmp/5mcic1321793415.ps tmp/5mcic1321793415.png",intern=TRUE))
character(0)
> try(system("convert tmp/6c7na1321793415.ps tmp/6c7na1321793415.png",intern=TRUE))
character(0)
> try(system("convert tmp/7tz8a1321793415.ps tmp/7tz8a1321793415.png",intern=TRUE))
character(0)
> try(system("convert tmp/8ch1d1321793415.ps tmp/8ch1d1321793415.png",intern=TRUE))
character(0)
> try(system("convert tmp/9ihdq1321793415.ps tmp/9ihdq1321793415.png",intern=TRUE))
character(0)
> try(system("convert tmp/10c7vl1321793415.ps tmp/10c7vl1321793415.png",intern=TRUE))
character(0)
>
>
> proc.time()
user system elapsed
3.190 0.445 3.743