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(300
+ ,2.26
+ ,591.000
+ ,302
+ ,2.57
+ ,589.000
+ ,400
+ ,3.07
+ ,584.000
+ ,392
+ ,2.76
+ ,573.000
+ ,373
+ ,2.51
+ ,567.000
+ ,379
+ ,2.87
+ ,569.000
+ ,303
+ ,3.14
+ ,621.000
+ ,324
+ ,3.11
+ ,629.000
+ ,353
+ ,3.16
+ ,628.000
+ ,392
+ ,2.47
+ ,612.000
+ ,327
+ ,2.57
+ ,595.000
+ ,376
+ ,2.89
+ ,597.000
+ ,329
+ ,2.63
+ ,593.000
+ ,359
+ ,2.38
+ ,590.000
+ ,413
+ ,1.69
+ ,580.000
+ ,338
+ ,1.96
+ ,574.000
+ ,422
+ ,2.19
+ ,573.000
+ ,390
+ ,1.87
+ ,573.000
+ ,370
+ ,1.60
+ ,620.000
+ ,367
+ ,1.63
+ ,626.000
+ ,406
+ ,1.22
+ ,620.000
+ ,418
+ ,1.21
+ ,588.000
+ ,346
+ ,1.49
+ ,566.000
+ ,350
+ ,1.64
+ ,557.000
+ ,330
+ ,1.66
+ ,561.000
+ ,318
+ ,1.77
+ ,549.000
+ ,382
+ ,1.82
+ ,532.000
+ ,337
+ ,1.78
+ ,526.000
+ ,372
+ ,1.28
+ ,511.000
+ ,422
+ ,1.29
+ ,499.000
+ ,428
+ ,1.37
+ ,555.000
+ ,426
+ ,1.12
+ ,565.000
+ ,396
+ ,1.51
+ ,542.000
+ ,458
+ ,2.24
+ ,527.000
+ ,315
+ ,2.94
+ ,510.000
+ ,337
+ ,3.09
+ ,514.000
+ ,386
+ ,3.46
+ ,517.000
+ ,352
+ ,3.64
+ ,508.000
+ ,383
+ ,4.39
+ ,493.000
+ ,439
+ ,4.15
+ ,490.000
+ ,397
+ ,5.21
+ ,469.000
+ ,453
+ ,5.80
+ ,478.000
+ ,363
+ ,5.91
+ ,528.000
+ ,365
+ ,5.39
+ ,534.000
+ ,474
+ ,5.46
+ ,518.000
+ ,373
+ ,4.72
+ ,506.000
+ ,403
+ ,3.14
+ ,502.000
+ ,384
+ ,2.63
+ ,516.000
+ ,364
+ ,2.32
+ ,528.000
+ ,361
+ ,1.93
+ ,533.000
+ ,419
+ ,0.62
+ ,536.000
+ ,352
+ ,0.60
+ ,537.000
+ ,363
+ ,-0.37
+ ,524.000
+ ,410
+ ,-1.10
+ ,536.000
+ ,361
+ ,-1.68
+ ,587.000
+ ,383
+ ,-0.78
+ ,597.000
+ ,342
+ ,-1.19
+ ,581.000
+ ,369
+ ,-0.79
+ ,564.000
+ ,361
+ ,-0.12
+ ,558.000
+ ,317
+ ,0.26
+ ,575.000
+ ,386
+ ,0.62
+ ,580.000
+ ,318
+ ,0.70
+ ,575.000
+ ,407
+ ,1.66
+ ,563.000
+ ,393
+ ,1.80
+ ,552.000
+ ,404
+ ,2.27
+ ,537.000
+ ,498
+ ,2.46
+ ,545.000
+ ,438
+ ,2.57
+ ,601.000)
+ ,dim=c(3
+ ,67)
+ ,dimnames=list(c('Aantal_vergunningen'
+ ,'Inflatie'
+ ,'Aantal_werklozen')
+ ,1:67))
> y <- array(NA,dim=c(3,67),dimnames=list(c('Aantal_vergunningen','Inflatie','Aantal_werklozen'),1:67))
> 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 = '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
Aantal_vergunningen Inflatie Aantal_werklozen
1 300 2.26 591
2 302 2.57 589
3 400 3.07 584
4 392 2.76 573
5 373 2.51 567
6 379 2.87 569
7 303 3.14 621
8 324 3.11 629
9 353 3.16 628
10 392 2.47 612
11 327 2.57 595
12 376 2.89 597
13 329 2.63 593
14 359 2.38 590
15 413 1.69 580
16 338 1.96 574
17 422 2.19 573
18 390 1.87 573
19 370 1.60 620
20 367 1.63 626
21 406 1.22 620
22 418 1.21 588
23 346 1.49 566
24 350 1.64 557
25 330 1.66 561
26 318 1.77 549
27 382 1.82 532
28 337 1.78 526
29 372 1.28 511
30 422 1.29 499
31 428 1.37 555
32 426 1.12 565
33 396 1.51 542
34 458 2.24 527
35 315 2.94 510
36 337 3.09 514
37 386 3.46 517
38 352 3.64 508
39 383 4.39 493
40 439 4.15 490
41 397 5.21 469
42 453 5.80 478
43 363 5.91 528
44 365 5.39 534
45 474 5.46 518
46 373 4.72 506
47 403 3.14 502
48 384 2.63 516
49 364 2.32 528
50 361 1.93 533
51 419 0.62 536
52 352 0.60 537
53 363 -0.37 524
54 410 -1.10 536
55 361 -1.68 587
56 383 -0.78 597
57 342 -1.19 581
58 369 -0.79 564
59 361 -0.12 558
60 317 0.26 575
61 386 0.62 580
62 318 0.70 575
63 407 1.66 563
64 393 1.80 552
65 404 2.27 537
66 498 2.46 545
67 438 2.57 601
> k <- length(x[1,])
> df <- as.data.frame(x)
> (mylm <- lm(df))
Call:
lm(formula = df)
Coefficients:
(Intercept) Inflatie Aantal_werklozen
540.4959 1.2806 -0.2985
> (mysum <- summary(mylm))
Call:
lm(formula = df)
Residuals:
Min 1Q Median 3Q Max
-77.01 -27.43 -2.11 29.37 117.06
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 540.4959 79.3694 6.810 4.01e-09 ***
Inflatie 1.2806 3.3226 0.385 0.7012
Aantal_werklozen -0.2985 0.1374 -2.173 0.0335 *
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 40.86 on 64 degrees of freedom
Multiple R-squared: 0.08883, Adjusted R-squared: 0.06036
F-statistic: 3.12 on 2 and 64 DF, p-value: 0.05095
> 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.16477033 0.3295407 0.8352297
[2,] 0.07538656 0.1507731 0.9246134
[3,] 0.07646861 0.1529372 0.9235314
[4,] 0.08949543 0.1789909 0.9105046
[5,] 0.53038138 0.9392372 0.4696186
[6,] 0.45162776 0.9032555 0.5483722
[7,] 0.37039564 0.7407913 0.6296044
[8,] 0.31542405 0.6308481 0.6845760
[9,] 0.25768537 0.5153707 0.7423146
[10,] 0.38655471 0.7731094 0.6134453
[11,] 0.36181846 0.7236369 0.6381815
[12,] 0.40653535 0.8130707 0.5934647
[13,] 0.33275735 0.6655147 0.6672426
[14,] 0.28067995 0.5613599 0.7193200
[15,] 0.22385732 0.4477146 0.7761427
[16,] 0.20239201 0.4047840 0.7976080
[17,] 0.17294304 0.3458861 0.8270570
[18,] 0.19915808 0.3983162 0.8008419
[19,] 0.19181402 0.3836280 0.8081860
[20,] 0.22802603 0.4560521 0.7719740
[21,] 0.29831925 0.5966385 0.7016808
[22,] 0.24143379 0.4828676 0.7585662
[23,] 0.23387949 0.4677590 0.7661205
[24,] 0.18184837 0.3636967 0.8181516
[25,] 0.19488457 0.3897691 0.8051154
[26,] 0.21400305 0.4280061 0.7859969
[27,] 0.21472545 0.4294509 0.7852745
[28,] 0.17018436 0.3403687 0.8298156
[29,] 0.34069966 0.6813993 0.6593003
[30,] 0.43283295 0.8656659 0.5671671
[31,] 0.43871825 0.8774365 0.5612818
[32,] 0.40707644 0.8141529 0.5929236
[33,] 0.39360809 0.7872162 0.6063919
[34,] 0.37329552 0.7465910 0.6267045
[35,] 0.43797927 0.8759585 0.5620207
[36,] 0.38651472 0.7730294 0.6134853
[37,] 0.44134518 0.8826904 0.5586548
[38,] 0.45277082 0.9055416 0.5472292
[39,] 0.50395730 0.9920854 0.4960427
[40,] 0.61538784 0.7692243 0.3846122
[41,] 0.61789226 0.7642155 0.3821077
[42,] 0.54009335 0.9198133 0.4599067
[43,] 0.47613296 0.9522659 0.5238670
[44,] 0.47637265 0.9527453 0.5236274
[45,] 0.50377171 0.9924566 0.4962283
[46,] 0.46310221 0.9262044 0.5368978
[47,] 0.45270107 0.9054021 0.5472989
[48,] 0.38781304 0.7756261 0.6121870
[49,] 0.41605064 0.8321013 0.5839494
[50,] 0.37625010 0.7525002 0.6237499
[51,] 0.39839643 0.7967929 0.6016036
[52,] 0.34372630 0.6874526 0.6562737
[53,] 0.42379697 0.8475939 0.5762030
[54,] 0.44151123 0.8830225 0.5584888
[55,] 0.32545360 0.6509072 0.6745464
[56,] 0.44392295 0.8878459 0.5560770
> postscript(file="/var/www/html/rcomp/tmp/1uyl71292940626.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/html/rcomp/tmp/2npka1292940626.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/html/rcomp/tmp/3npka1292940626.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/html/rcomp/tmp/4npka1292940626.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/html/rcomp/tmp/5ghkd1292940626.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 = 67
Frequency = 1
1 2 3 4 5 6 7
-66.955234 -65.949282 29.917754 19.030836 -1.440237 4.695830 -56.126059
8 9 10 11 12 13 14
-32.699354 -4.061919 31.045100 -39.158067 10.029223 -37.831973 -8.407438
15 16 17 18 19 20 21
43.490796 -33.646173 49.760760 18.170542 12.547482 11.300280 49.034098
22 23 24 25 26 27 28
51.493756 -27.432593 -26.311501 -45.142969 -60.866262 -2.005400 -48.745393
29 30 31 32 33 34 35
-17.583146 28.821618 51.437181 52.742682 15.376935 71.964081 -77.007427
36 37 38 39 40 41 42
-54.005369 -4.583572 -41.500897 -15.939362 39.472367 -10.154290 47.776997
43 44 45 46 47 48 49
-27.437072 -22.979960 81.153826 -22.480984 8.348172 -5.819235 -21.839828
50 51 52 53 54 55 56
-22.847727 37.725426 -28.950426 -20.589241 30.928005 -2.103935 21.728912
57 58 59 60 61 62 63
-23.522629 -2.109966 -12.759163 -52.170669 17.861005 -51.734120 32.454103
64 65 66 67
14.990929 20.911023 117.056002 73.633148
> postscript(file="/var/www/html/rcomp/tmp/6ghkd1292940626.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 = 67
Frequency = 1
lag(myerror, k = 1) myerror
0 -66.955234 NA
1 -65.949282 -66.955234
2 29.917754 -65.949282
3 19.030836 29.917754
4 -1.440237 19.030836
5 4.695830 -1.440237
6 -56.126059 4.695830
7 -32.699354 -56.126059
8 -4.061919 -32.699354
9 31.045100 -4.061919
10 -39.158067 31.045100
11 10.029223 -39.158067
12 -37.831973 10.029223
13 -8.407438 -37.831973
14 43.490796 -8.407438
15 -33.646173 43.490796
16 49.760760 -33.646173
17 18.170542 49.760760
18 12.547482 18.170542
19 11.300280 12.547482
20 49.034098 11.300280
21 51.493756 49.034098
22 -27.432593 51.493756
23 -26.311501 -27.432593
24 -45.142969 -26.311501
25 -60.866262 -45.142969
26 -2.005400 -60.866262
27 -48.745393 -2.005400
28 -17.583146 -48.745393
29 28.821618 -17.583146
30 51.437181 28.821618
31 52.742682 51.437181
32 15.376935 52.742682
33 71.964081 15.376935
34 -77.007427 71.964081
35 -54.005369 -77.007427
36 -4.583572 -54.005369
37 -41.500897 -4.583572
38 -15.939362 -41.500897
39 39.472367 -15.939362
40 -10.154290 39.472367
41 47.776997 -10.154290
42 -27.437072 47.776997
43 -22.979960 -27.437072
44 81.153826 -22.979960
45 -22.480984 81.153826
46 8.348172 -22.480984
47 -5.819235 8.348172
48 -21.839828 -5.819235
49 -22.847727 -21.839828
50 37.725426 -22.847727
51 -28.950426 37.725426
52 -20.589241 -28.950426
53 30.928005 -20.589241
54 -2.103935 30.928005
55 21.728912 -2.103935
56 -23.522629 21.728912
57 -2.109966 -23.522629
58 -12.759163 -2.109966
59 -52.170669 -12.759163
60 17.861005 -52.170669
61 -51.734120 17.861005
62 32.454103 -51.734120
63 14.990929 32.454103
64 20.911023 14.990929
65 117.056002 20.911023
66 73.633148 117.056002
67 NA 73.633148
> dum1 <- dum[2:length(myerror),]
> dum1
lag(myerror, k = 1) myerror
[1,] -65.949282 -66.955234
[2,] 29.917754 -65.949282
[3,] 19.030836 29.917754
[4,] -1.440237 19.030836
[5,] 4.695830 -1.440237
[6,] -56.126059 4.695830
[7,] -32.699354 -56.126059
[8,] -4.061919 -32.699354
[9,] 31.045100 -4.061919
[10,] -39.158067 31.045100
[11,] 10.029223 -39.158067
[12,] -37.831973 10.029223
[13,] -8.407438 -37.831973
[14,] 43.490796 -8.407438
[15,] -33.646173 43.490796
[16,] 49.760760 -33.646173
[17,] 18.170542 49.760760
[18,] 12.547482 18.170542
[19,] 11.300280 12.547482
[20,] 49.034098 11.300280
[21,] 51.493756 49.034098
[22,] -27.432593 51.493756
[23,] -26.311501 -27.432593
[24,] -45.142969 -26.311501
[25,] -60.866262 -45.142969
[26,] -2.005400 -60.866262
[27,] -48.745393 -2.005400
[28,] -17.583146 -48.745393
[29,] 28.821618 -17.583146
[30,] 51.437181 28.821618
[31,] 52.742682 51.437181
[32,] 15.376935 52.742682
[33,] 71.964081 15.376935
[34,] -77.007427 71.964081
[35,] -54.005369 -77.007427
[36,] -4.583572 -54.005369
[37,] -41.500897 -4.583572
[38,] -15.939362 -41.500897
[39,] 39.472367 -15.939362
[40,] -10.154290 39.472367
[41,] 47.776997 -10.154290
[42,] -27.437072 47.776997
[43,] -22.979960 -27.437072
[44,] 81.153826 -22.979960
[45,] -22.480984 81.153826
[46,] 8.348172 -22.480984
[47,] -5.819235 8.348172
[48,] -21.839828 -5.819235
[49,] -22.847727 -21.839828
[50,] 37.725426 -22.847727
[51,] -28.950426 37.725426
[52,] -20.589241 -28.950426
[53,] 30.928005 -20.589241
[54,] -2.103935 30.928005
[55,] 21.728912 -2.103935
[56,] -23.522629 21.728912
[57,] -2.109966 -23.522629
[58,] -12.759163 -2.109966
[59,] -52.170669 -12.759163
[60,] 17.861005 -52.170669
[61,] -51.734120 17.861005
[62,] 32.454103 -51.734120
[63,] 14.990929 32.454103
[64,] 20.911023 14.990929
[65,] 117.056002 20.911023
[66,] 73.633148 117.056002
> z <- as.data.frame(dum1)
> z
lag(myerror, k = 1) myerror
1 -65.949282 -66.955234
2 29.917754 -65.949282
3 19.030836 29.917754
4 -1.440237 19.030836
5 4.695830 -1.440237
6 -56.126059 4.695830
7 -32.699354 -56.126059
8 -4.061919 -32.699354
9 31.045100 -4.061919
10 -39.158067 31.045100
11 10.029223 -39.158067
12 -37.831973 10.029223
13 -8.407438 -37.831973
14 43.490796 -8.407438
15 -33.646173 43.490796
16 49.760760 -33.646173
17 18.170542 49.760760
18 12.547482 18.170542
19 11.300280 12.547482
20 49.034098 11.300280
21 51.493756 49.034098
22 -27.432593 51.493756
23 -26.311501 -27.432593
24 -45.142969 -26.311501
25 -60.866262 -45.142969
26 -2.005400 -60.866262
27 -48.745393 -2.005400
28 -17.583146 -48.745393
29 28.821618 -17.583146
30 51.437181 28.821618
31 52.742682 51.437181
32 15.376935 52.742682
33 71.964081 15.376935
34 -77.007427 71.964081
35 -54.005369 -77.007427
36 -4.583572 -54.005369
37 -41.500897 -4.583572
38 -15.939362 -41.500897
39 39.472367 -15.939362
40 -10.154290 39.472367
41 47.776997 -10.154290
42 -27.437072 47.776997
43 -22.979960 -27.437072
44 81.153826 -22.979960
45 -22.480984 81.153826
46 8.348172 -22.480984
47 -5.819235 8.348172
48 -21.839828 -5.819235
49 -22.847727 -21.839828
50 37.725426 -22.847727
51 -28.950426 37.725426
52 -20.589241 -28.950426
53 30.928005 -20.589241
54 -2.103935 30.928005
55 21.728912 -2.103935
56 -23.522629 21.728912
57 -2.109966 -23.522629
58 -12.759163 -2.109966
59 -52.170669 -12.759163
60 17.861005 -52.170669
61 -51.734120 17.861005
62 32.454103 -51.734120
63 14.990929 32.454103
64 20.911023 14.990929
65 117.056002 20.911023
66 73.633148 117.056002
> 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/7r81g1292940626.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/html/rcomp/tmp/8r81g1292940626.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/html/rcomp/tmp/91z0i1292940626.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/html/rcomp/tmp/101z0i1292940626.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/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/11n0zo1292940626.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/12q0fu1292940626.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/13xju61292940626.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/14psc91292940626.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/15tbax1292940626.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/16etr31292940626.tab")
+ }
>
> try(system("convert tmp/1uyl71292940626.ps tmp/1uyl71292940626.png",intern=TRUE))
character(0)
> try(system("convert tmp/2npka1292940626.ps tmp/2npka1292940626.png",intern=TRUE))
character(0)
> try(system("convert tmp/3npka1292940626.ps tmp/3npka1292940626.png",intern=TRUE))
character(0)
> try(system("convert tmp/4npka1292940626.ps tmp/4npka1292940626.png",intern=TRUE))
character(0)
> try(system("convert tmp/5ghkd1292940626.ps tmp/5ghkd1292940626.png",intern=TRUE))
character(0)
> try(system("convert tmp/6ghkd1292940626.ps tmp/6ghkd1292940626.png",intern=TRUE))
character(0)
> try(system("convert tmp/7r81g1292940626.ps tmp/7r81g1292940626.png",intern=TRUE))
character(0)
> try(system("convert tmp/8r81g1292940626.ps tmp/8r81g1292940626.png",intern=TRUE))
character(0)
> try(system("convert tmp/91z0i1292940626.ps tmp/91z0i1292940626.png",intern=TRUE))
character(0)
> try(system("convert tmp/101z0i1292940626.ps tmp/101z0i1292940626.png",intern=TRUE))
character(0)
>
>
> proc.time()
user system elapsed
2.561 1.620 5.815