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(631.923
+ ,9.911
+ ,58608
+ ,654.294
+ ,8.915
+ ,46865
+ ,671.833
+ ,9.452
+ ,51378
+ ,586.840
+ ,9.112
+ ,46235
+ ,600.969
+ ,8.472
+ ,47206
+ ,625.568
+ ,8.230
+ ,45382
+ ,558.110
+ ,8.384
+ ,41227
+ ,630.577
+ ,8.625
+ ,33795
+ ,628.654
+ ,8.221
+ ,31295
+ ,603.184
+ ,8.649
+ ,42625
+ ,656.255
+ ,8.625
+ ,33625
+ ,600.730
+ ,10.443
+ ,21538
+ ,670.326
+ ,10.357
+ ,56421
+ ,678.423
+ ,8.586
+ ,53152
+ ,641.502
+ ,8.892
+ ,53536
+ ,625.311
+ ,8.329
+ ,52408
+ ,628.177
+ ,8.101
+ ,41454
+ ,589.767
+ ,7.922
+ ,38271
+ ,582.471
+ ,8.120
+ ,35306
+ ,636.248
+ ,7.838
+ ,26414
+ ,599.885
+ ,7.735
+ ,31917
+ ,621.694
+ ,8.406
+ ,38030
+ ,637.406
+ ,8.209
+ ,27534
+ ,595.994
+ ,9.451
+ ,18387
+ ,696.308
+ ,10.041
+ ,50556
+ ,674.201
+ ,9.411
+ ,43901
+ ,648.861
+ ,10.405
+ ,48572
+ ,649.605
+ ,8.467
+ ,43899
+ ,672.392
+ ,8.464
+ ,37532
+ ,598.396
+ ,8.102
+ ,40357
+ ,613.177
+ ,7.627
+ ,35489
+ ,638.104
+ ,7.513
+ ,29027
+ ,615.632
+ ,7.510
+ ,34485
+ ,634.465
+ ,8.291
+ ,42598
+ ,638.686
+ ,8.064
+ ,30306
+ ,604.243
+ ,9.383
+ ,26451
+ ,706.669
+ ,9.706
+ ,47460
+ ,677.185
+ ,8.579
+ ,50104
+ ,644.328
+ ,9.474
+ ,61465
+ ,664.825
+ ,8.318
+ ,53726
+ ,605.707
+ ,8.213
+ ,39477
+ ,600.136
+ ,8.059
+ ,43895
+ ,612.166
+ ,9.111
+ ,31481
+ ,599.659
+ ,7.708
+ ,29896
+ ,634.210
+ ,7.680
+ ,33842
+ ,618.234
+ ,8.014
+ ,39120
+ ,613.576
+ ,8.007
+ ,33702
+ ,627.200
+ ,8.718
+ ,25094
+ ,668.973
+ ,9.486
+ ,51442
+ ,651.479
+ ,9.113
+ ,45594
+ ,619.661
+ ,9.025
+ ,52518
+ ,644.260
+ ,8.476
+ ,48564
+ ,579.936
+ ,7.952
+ ,41745
+ ,601.752
+ ,7.759
+ ,49585
+ ,595.376
+ ,7.835
+ ,32747
+ ,588.902
+ ,7.600
+ ,33379
+ ,634.341
+ ,7.651
+ ,35645
+ ,594.305
+ ,8.319
+ ,37034
+ ,606.200
+ ,8.812
+ ,35681
+ ,610.926
+ ,8.630
+ ,20972)
+ ,dim=c(3
+ ,60)
+ ,dimnames=list(c('WERKLOZEN'
+ ,'OVERLIJDENS'
+ ,'INSCHRIJVINGEN')
+ ,1:60))
> y <- array(NA,dim=c(3,60),dimnames=list(c('WERKLOZEN','OVERLIJDENS','INSCHRIJVINGEN'),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 = 'No Linear Trend'
> par2 = 'Do not include Seasonal Dummies'
> par1 = '2'
> 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
OVERLIJDENS WERKLOZEN INSCHRIJVINGEN
1 9.911 631.923 58608
2 8.915 654.294 46865
3 9.452 671.833 51378
4 9.112 586.840 46235
5 8.472 600.969 47206
6 8.230 625.568 45382
7 8.384 558.110 41227
8 8.625 630.577 33795
9 8.221 628.654 31295
10 8.649 603.184 42625
11 8.625 656.255 33625
12 10.443 600.730 21538
13 10.357 670.326 56421
14 8.586 678.423 53152
15 8.892 641.502 53536
16 8.329 625.311 52408
17 8.101 628.177 41454
18 7.922 589.767 38271
19 8.120 582.471 35306
20 7.838 636.248 26414
21 7.735 599.885 31917
22 8.406 621.694 38030
23 8.209 637.406 27534
24 9.451 595.994 18387
25 10.041 696.308 50556
26 9.411 674.201 43901
27 10.405 648.861 48572
28 8.467 649.605 43899
29 8.464 672.392 37532
30 8.102 598.396 40357
31 7.627 613.177 35489
32 7.513 638.104 29027
33 7.510 615.632 34485
34 8.291 634.465 42598
35 8.064 638.686 30306
36 9.383 604.243 26451
37 9.706 706.669 47460
38 8.579 677.185 50104
39 9.474 644.328 61465
40 8.318 664.825 53726
41 8.213 605.707 39477
42 8.059 600.136 43895
43 9.111 612.166 31481
44 7.708 599.659 29896
45 7.680 634.210 33842
46 8.014 618.234 39120
47 8.007 613.576 33702
48 8.718 627.200 25094
49 9.486 668.973 51442
50 9.113 651.479 45594
51 9.025 619.661 52518
52 8.476 644.260 48564
53 7.952 579.936 41745
54 7.759 601.752 49585
55 7.835 595.376 32747
56 7.600 588.902 33379
57 7.651 634.341 35645
58 8.319 594.305 37034
59 8.812 606.200 35681
60 8.630 610.926 20972
> k <- length(x[1,])
> df <- as.data.frame(x)
> (mylm <- lm(df))
Call:
lm(formula = df)
Coefficients:
(Intercept) WERKLOZEN INSCHRIJVINGEN
2.870e+00 8.157e-03 1.463e-05
> (mysum <- summary(mylm))
Call:
lm(formula = df)
Residuals:
Min 1Q Median 3Q Max
-0.9870 -0.4625 -0.2053 0.3676 2.3575
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 2.870e+00 1.862e+00 1.541 0.1288
WERKLOZEN 8.157e-03 3.185e-03 2.561 0.0131 *
INSCHRIJVINGEN 1.463e-05 9.980e-06 1.466 0.1481
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 0.6848 on 57 degrees of freedom
Multiple R-squared: 0.2053, Adjusted R-squared: 0.1774
F-statistic: 7.361 on 2 and 57 DF, p-value: 0.001434
> 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.20414619 0.40829238 0.79585381
[2,] 0.10547069 0.21094138 0.89452931
[3,] 0.13977430 0.27954860 0.86022570
[4,] 0.07276225 0.14552449 0.92723775
[5,] 0.03521988 0.07043975 0.96478012
[6,] 0.01768368 0.03536736 0.98231632
[7,] 0.80040433 0.39919134 0.19959567
[8,] 0.89446968 0.21106065 0.10553032
[9,] 0.89144213 0.21711574 0.10855787
[10,] 0.84579266 0.30841467 0.15420734
[11,] 0.82420634 0.35158733 0.17579366
[12,] 0.82452569 0.35094862 0.17547431
[13,] 0.81008266 0.37983468 0.18991734
[14,] 0.76353308 0.47293384 0.23646692
[15,] 0.77491037 0.45017925 0.22508963
[16,] 0.75711446 0.48577107 0.24288554
[17,] 0.69288799 0.61422403 0.30711201
[18,] 0.63121107 0.73757786 0.36878893
[19,] 0.82137689 0.35724621 0.17862311
[20,] 0.82924447 0.34151105 0.17075553
[21,] 0.79215831 0.41568338 0.20784169
[22,] 0.94978304 0.10043392 0.05021696
[23,] 0.93425278 0.13149444 0.06574722
[24,] 0.91901216 0.16197568 0.08098784
[25,] 0.89264377 0.21471246 0.10735623
[26,] 0.90073707 0.19852585 0.09926293
[27,] 0.93419245 0.13161511 0.06580755
[28,] 0.94915011 0.10169978 0.05084989
[29,] 0.93147580 0.13704840 0.06852420
[30,] 0.91890237 0.16219526 0.08109763
[31,] 0.97163525 0.05672950 0.02836475
[32,] 0.96187690 0.07624621 0.03812310
[33,] 0.95421873 0.09156254 0.04578127
[34,] 0.95602566 0.08794867 0.04397434
[35,] 0.95864609 0.08270782 0.04135391
[36,] 0.93564739 0.12870521 0.06435261
[37,] 0.90583247 0.18833506 0.09416753
[38,] 0.93963516 0.12072969 0.06036484
[39,] 0.91673728 0.16652544 0.08326272
[40,] 0.94881244 0.10237513 0.05118756
[41,] 0.93115711 0.13768577 0.06884289
[42,] 0.90221134 0.19557733 0.09778866
[43,] 0.85390945 0.29218110 0.14609055
[44,] 0.80601839 0.38796321 0.19398161
[45,] 0.75418208 0.49163583 0.24581792
[46,] 0.82241978 0.35516044 0.17758022
[47,] 0.78058994 0.43882012 0.21941006
[48,] 0.64730709 0.70538582 0.35269291
[49,] 0.48578197 0.97156395 0.51421803
> postscript(file="/var/www/rcomp/tmp/1qw2o1324637792.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/2qrfk1324637792.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/3fga51324637792.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/4qxi41324637792.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/59m9e1324637792.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
1.028608067 0.021943199 0.349835140 0.778413066 0.008948867 -0.407026871
7 8 9 10 11 12
0.358053886 0.116655036 -0.235077890 0.234909939 -0.090324119 2.357476027
13 14 15 16 17 18
1.193338559 -0.595879757 0.005682164 -0.408735762 -0.499834610 -0.318933457
19 20 21 22 23 24
-0.018032543 -0.608606021 -0.495497826 -0.091849543 -0.263440317 1.450215483
25 26 27 28 29 30
0.751209528 0.398922775 1.531285497 -0.344407688 -0.440128384 -0.239846589
31 32 33 34 35 36
-0.764192424 -0.986979978 -0.886528268 -0.377867719 -0.459442101 1.196931299
37 38 39 40 41 42
0.376991421 -0.548182118 0.448611087 -0.761353799 -0.175609322 -0.348809005
43 44 45 46 47 48
0.786700323 -0.491082709 -0.858668636 -0.471573840 -0.361299623 0.364516864
49 50 51 52 53 54
0.406228948 0.261503817 0.331744212 -0.360065088 -0.259569736 -0.745248244
55 56 57 58 59 60
-0.370860620 -0.562296892 -0.915118992 0.059148048 0.474912631 0.469584605
> postscript(file="/var/www/rcomp/tmp/6wqj01324637792.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 1.028608067 NA
1 0.021943199 1.028608067
2 0.349835140 0.021943199
3 0.778413066 0.349835140
4 0.008948867 0.778413066
5 -0.407026871 0.008948867
6 0.358053886 -0.407026871
7 0.116655036 0.358053886
8 -0.235077890 0.116655036
9 0.234909939 -0.235077890
10 -0.090324119 0.234909939
11 2.357476027 -0.090324119
12 1.193338559 2.357476027
13 -0.595879757 1.193338559
14 0.005682164 -0.595879757
15 -0.408735762 0.005682164
16 -0.499834610 -0.408735762
17 -0.318933457 -0.499834610
18 -0.018032543 -0.318933457
19 -0.608606021 -0.018032543
20 -0.495497826 -0.608606021
21 -0.091849543 -0.495497826
22 -0.263440317 -0.091849543
23 1.450215483 -0.263440317
24 0.751209528 1.450215483
25 0.398922775 0.751209528
26 1.531285497 0.398922775
27 -0.344407688 1.531285497
28 -0.440128384 -0.344407688
29 -0.239846589 -0.440128384
30 -0.764192424 -0.239846589
31 -0.986979978 -0.764192424
32 -0.886528268 -0.986979978
33 -0.377867719 -0.886528268
34 -0.459442101 -0.377867719
35 1.196931299 -0.459442101
36 0.376991421 1.196931299
37 -0.548182118 0.376991421
38 0.448611087 -0.548182118
39 -0.761353799 0.448611087
40 -0.175609322 -0.761353799
41 -0.348809005 -0.175609322
42 0.786700323 -0.348809005
43 -0.491082709 0.786700323
44 -0.858668636 -0.491082709
45 -0.471573840 -0.858668636
46 -0.361299623 -0.471573840
47 0.364516864 -0.361299623
48 0.406228948 0.364516864
49 0.261503817 0.406228948
50 0.331744212 0.261503817
51 -0.360065088 0.331744212
52 -0.259569736 -0.360065088
53 -0.745248244 -0.259569736
54 -0.370860620 -0.745248244
55 -0.562296892 -0.370860620
56 -0.915118992 -0.562296892
57 0.059148048 -0.915118992
58 0.474912631 0.059148048
59 0.469584605 0.474912631
60 NA 0.469584605
> dum1 <- dum[2:length(myerror),]
> dum1
lag(myerror, k = 1) myerror
[1,] 0.021943199 1.028608067
[2,] 0.349835140 0.021943199
[3,] 0.778413066 0.349835140
[4,] 0.008948867 0.778413066
[5,] -0.407026871 0.008948867
[6,] 0.358053886 -0.407026871
[7,] 0.116655036 0.358053886
[8,] -0.235077890 0.116655036
[9,] 0.234909939 -0.235077890
[10,] -0.090324119 0.234909939
[11,] 2.357476027 -0.090324119
[12,] 1.193338559 2.357476027
[13,] -0.595879757 1.193338559
[14,] 0.005682164 -0.595879757
[15,] -0.408735762 0.005682164
[16,] -0.499834610 -0.408735762
[17,] -0.318933457 -0.499834610
[18,] -0.018032543 -0.318933457
[19,] -0.608606021 -0.018032543
[20,] -0.495497826 -0.608606021
[21,] -0.091849543 -0.495497826
[22,] -0.263440317 -0.091849543
[23,] 1.450215483 -0.263440317
[24,] 0.751209528 1.450215483
[25,] 0.398922775 0.751209528
[26,] 1.531285497 0.398922775
[27,] -0.344407688 1.531285497
[28,] -0.440128384 -0.344407688
[29,] -0.239846589 -0.440128384
[30,] -0.764192424 -0.239846589
[31,] -0.986979978 -0.764192424
[32,] -0.886528268 -0.986979978
[33,] -0.377867719 -0.886528268
[34,] -0.459442101 -0.377867719
[35,] 1.196931299 -0.459442101
[36,] 0.376991421 1.196931299
[37,] -0.548182118 0.376991421
[38,] 0.448611087 -0.548182118
[39,] -0.761353799 0.448611087
[40,] -0.175609322 -0.761353799
[41,] -0.348809005 -0.175609322
[42,] 0.786700323 -0.348809005
[43,] -0.491082709 0.786700323
[44,] -0.858668636 -0.491082709
[45,] -0.471573840 -0.858668636
[46,] -0.361299623 -0.471573840
[47,] 0.364516864 -0.361299623
[48,] 0.406228948 0.364516864
[49,] 0.261503817 0.406228948
[50,] 0.331744212 0.261503817
[51,] -0.360065088 0.331744212
[52,] -0.259569736 -0.360065088
[53,] -0.745248244 -0.259569736
[54,] -0.370860620 -0.745248244
[55,] -0.562296892 -0.370860620
[56,] -0.915118992 -0.562296892
[57,] 0.059148048 -0.915118992
[58,] 0.474912631 0.059148048
[59,] 0.469584605 0.474912631
> z <- as.data.frame(dum1)
> z
lag(myerror, k = 1) myerror
1 0.021943199 1.028608067
2 0.349835140 0.021943199
3 0.778413066 0.349835140
4 0.008948867 0.778413066
5 -0.407026871 0.008948867
6 0.358053886 -0.407026871
7 0.116655036 0.358053886
8 -0.235077890 0.116655036
9 0.234909939 -0.235077890
10 -0.090324119 0.234909939
11 2.357476027 -0.090324119
12 1.193338559 2.357476027
13 -0.595879757 1.193338559
14 0.005682164 -0.595879757
15 -0.408735762 0.005682164
16 -0.499834610 -0.408735762
17 -0.318933457 -0.499834610
18 -0.018032543 -0.318933457
19 -0.608606021 -0.018032543
20 -0.495497826 -0.608606021
21 -0.091849543 -0.495497826
22 -0.263440317 -0.091849543
23 1.450215483 -0.263440317
24 0.751209528 1.450215483
25 0.398922775 0.751209528
26 1.531285497 0.398922775
27 -0.344407688 1.531285497
28 -0.440128384 -0.344407688
29 -0.239846589 -0.440128384
30 -0.764192424 -0.239846589
31 -0.986979978 -0.764192424
32 -0.886528268 -0.986979978
33 -0.377867719 -0.886528268
34 -0.459442101 -0.377867719
35 1.196931299 -0.459442101
36 0.376991421 1.196931299
37 -0.548182118 0.376991421
38 0.448611087 -0.548182118
39 -0.761353799 0.448611087
40 -0.175609322 -0.761353799
41 -0.348809005 -0.175609322
42 0.786700323 -0.348809005
43 -0.491082709 0.786700323
44 -0.858668636 -0.491082709
45 -0.471573840 -0.858668636
46 -0.361299623 -0.471573840
47 0.364516864 -0.361299623
48 0.406228948 0.364516864
49 0.261503817 0.406228948
50 0.331744212 0.261503817
51 -0.360065088 0.331744212
52 -0.259569736 -0.360065088
53 -0.745248244 -0.259569736
54 -0.370860620 -0.745248244
55 -0.562296892 -0.370860620
56 -0.915118992 -0.562296892
57 0.059148048 -0.915118992
58 0.474912631 0.059148048
59 0.469584605 0.474912631
> 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/7lpty1324637792.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/8e7341324637792.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/936tu1324637792.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/10x8fc1324637792.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/11pphf1324637792.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/12n0kx1324637792.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/13ntbj1324637792.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/14juk51324637792.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/15kd3y1324637792.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/16k73o1324637792.tab")
+ }
>
> try(system("convert tmp/1qw2o1324637792.ps tmp/1qw2o1324637792.png",intern=TRUE))
character(0)
> try(system("convert tmp/2qrfk1324637792.ps tmp/2qrfk1324637792.png",intern=TRUE))
character(0)
> try(system("convert tmp/3fga51324637792.ps tmp/3fga51324637792.png",intern=TRUE))
character(0)
> try(system("convert tmp/4qxi41324637792.ps tmp/4qxi41324637792.png",intern=TRUE))
character(0)
> try(system("convert tmp/59m9e1324637792.ps tmp/59m9e1324637792.png",intern=TRUE))
character(0)
> try(system("convert tmp/6wqj01324637792.ps tmp/6wqj01324637792.png",intern=TRUE))
character(0)
> try(system("convert tmp/7lpty1324637792.ps tmp/7lpty1324637792.png",intern=TRUE))
character(0)
> try(system("convert tmp/8e7341324637792.ps tmp/8e7341324637792.png",intern=TRUE))
character(0)
> try(system("convert tmp/936tu1324637792.ps tmp/936tu1324637792.png",intern=TRUE))
character(0)
> try(system("convert tmp/10x8fc1324637792.ps tmp/10x8fc1324637792.png",intern=TRUE))
character(0)
>
>
> proc.time()
user system elapsed
3.430 0.190 3.601