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(589
+ ,248.85
+ ,65453
+ ,559
+ ,249.68
+ ,65715
+ ,623
+ ,251.13
+ ,66261
+ ,617
+ ,251.24
+ ,66332
+ ,603
+ ,253.24
+ ,66229
+ ,558
+ ,254.66
+ ,66579
+ ,609
+ ,255.85
+ ,66817
+ ,583
+ ,256.93
+ ,67373
+ ,570
+ ,258.99
+ ,68078
+ ,543
+ ,258.30
+ ,69137
+ ,598
+ ,260.53
+ ,69816
+ ,569
+ ,260.65
+ ,70252
+ ,552
+ ,260.98
+ ,70389
+ ,514
+ ,262.09
+ ,70572
+ ,569
+ ,263.18
+ ,70780
+ ,529
+ ,262.62
+ ,70912
+ ,515
+ ,263.18
+ ,71594
+ ,481
+ ,264.91
+ ,72587
+ ,536
+ ,265.20
+ ,73677
+ ,498
+ ,266.14
+ ,74712
+ ,446
+ ,268.15
+ ,75208
+ ,503
+ ,270.62
+ ,75657
+ ,470
+ ,272.65
+ ,76011
+ ,458
+ ,274.50
+ ,76748
+ ,437
+ ,274.37
+ ,76537
+ ,502
+ ,277.85
+ ,76622
+ ,482
+ ,280.15
+ ,76404
+ ,474
+ ,280.67
+ ,76219
+ ,457
+ ,281.42
+ ,76875
+ ,522
+ ,283.23
+ ,77374
+ ,513
+ ,283.34
+ ,77743
+ ,515
+ ,284.09
+ ,78030
+ ,506
+ ,285.47
+ ,77805
+ ,576
+ ,287.27
+ ,77905
+ ,556
+ ,287.96
+ ,78158
+ ,559
+ ,289.05
+ ,78616
+ ,541
+ ,289.84
+ ,79740
+ ,606
+ ,292.68
+ ,80312
+ ,600
+ ,294.61
+ ,80921
+ ,588
+ ,296.22
+ ,81078
+ ,570
+ ,296.70
+ ,81394
+ ,626
+ ,300.82
+ ,81787
+ ,601
+ ,303.57
+ ,82252
+ ,588
+ ,304.32
+ ,82854
+ ,573
+ ,304.52
+ ,83498
+ ,622
+ ,306.69
+ ,83811
+ ,570
+ ,308.73
+ ,84531
+ ,547
+ ,308.30
+ ,85330
+ ,512
+ ,309.67
+ ,86247
+ ,554
+ ,311.68
+ ,86386
+ ,517
+ ,312.62
+ ,86918
+ ,506
+ ,315.18
+ ,87184
+ ,479
+ ,320.19
+ ,87843
+ ,527
+ ,325.96
+ ,88204
+ ,508
+ ,330.45
+ ,87675
+ ,532
+ ,329.16
+ ,85964
+ ,532
+ ,327.53
+ ,84387
+ ,588
+ ,326.87
+ ,84530
+ ,566
+ ,326.52
+ ,85497
+ ,573
+ ,326.65
+ ,85968
+ ,545
+ ,329.25
+ ,86030
+ ,597
+ ,333.11
+ ,86963
+ ,555
+ ,334.51
+ ,87324
+ ,548
+ ,336.21
+ ,87770
+ ,524
+ ,339.91
+ ,88534
+ ,572
+ ,344.53
+ ,88888)
+ ,dim=c(3
+ ,66)
+ ,dimnames=list(c('Werkloosheid'
+ ,'CPI'
+ ,'BBP')
+ ,1:66))
> y <- array(NA,dim=c(3,66),dimnames=list(c('Werkloosheid','CPI','BBP'),1:66))
> 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
> 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
Werkloosheid CPI BBP
1 589 248.85 65453
2 559 249.68 65715
3 623 251.13 66261
4 617 251.24 66332
5 603 253.24 66229
6 558 254.66 66579
7 609 255.85 66817
8 583 256.93 67373
9 570 258.99 68078
10 543 258.30 69137
11 598 260.53 69816
12 569 260.65 70252
13 552 260.98 70389
14 514 262.09 70572
15 569 263.18 70780
16 529 262.62 70912
17 515 263.18 71594
18 481 264.91 72587
19 536 265.20 73677
20 498 266.14 74712
21 446 268.15 75208
22 503 270.62 75657
23 470 272.65 76011
24 458 274.50 76748
25 437 274.37 76537
26 502 277.85 76622
27 482 280.15 76404
28 474 280.67 76219
29 457 281.42 76875
30 522 283.23 77374
31 513 283.34 77743
32 515 284.09 78030
33 506 285.47 77805
34 576 287.27 77905
35 556 287.96 78158
36 559 289.05 78616
37 541 289.84 79740
38 606 292.68 80312
39 600 294.61 80921
40 588 296.22 81078
41 570 296.70 81394
42 626 300.82 81787
43 601 303.57 82252
44 588 304.32 82854
45 573 304.52 83498
46 622 306.69 83811
47 570 308.73 84531
48 547 308.30 85330
49 512 309.67 86247
50 554 311.68 86386
51 517 312.62 86918
52 506 315.18 87184
53 479 320.19 87843
54 527 325.96 88204
55 508 330.45 87675
56 532 329.16 85964
57 532 327.53 84387
58 588 326.87 84530
59 566 326.52 85497
60 573 326.65 85968
61 545 329.25 86030
62 597 333.11 86963
63 555 334.51 87324
64 548 336.21 87770
65 524 339.91 88534
66 572 344.53 88888
> k <- length(x[1,])
> df <- as.data.frame(x)
> (mylm <- lm(df))
Call:
lm(formula = df)
Coefficients:
(Intercept) CPI BBP
592.38013 2.86224 -0.01121
> (mysum <- summary(mylm))
Call:
lm(formula = df)
Residuals:
Min 1Q Median 3Q Max
-82.540 -26.206 -1.353 28.791 91.511
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 592.380127 56.712918 10.445 2.23e-15 ***
CPI 2.862243 0.722471 3.962 0.000192 ***
BBP -0.011212 0.002775 -4.040 0.000148 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 42.14 on 63 degrees of freedom
Multiple R-squared: 0.2063, Adjusted R-squared: 0.1811
F-statistic: 8.188 on 2 and 63 DF, p-value: 0.0006902
> 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.282115156 0.56423031 0.71788484
[2,] 0.186748271 0.37349654 0.81325173
[3,] 0.155349553 0.31069911 0.84465045
[4,] 0.126404663 0.25280933 0.87359534
[5,] 0.097015192 0.19403038 0.90298481
[6,] 0.111021309 0.22204262 0.88897869
[7,] 0.069045583 0.13809117 0.93095442
[8,] 0.044778981 0.08955796 0.95522102
[9,] 0.058218270 0.11643654 0.94178173
[10,] 0.043885624 0.08777125 0.95611438
[11,] 0.031146131 0.06229226 0.96885387
[12,] 0.022003968 0.04400794 0.97799603
[13,] 0.022268694 0.04453739 0.97773131
[14,] 0.022890661 0.04578132 0.97710934
[15,] 0.013116478 0.02623296 0.98688352
[16,] 0.018073112 0.03614622 0.98192689
[17,] 0.012949658 0.02589932 0.98705034
[18,] 0.008974530 0.01794906 0.99102547
[19,] 0.007372019 0.01474404 0.99262798
[20,] 0.012914119 0.02582824 0.98708588
[21,] 0.012175041 0.02435008 0.98782496
[22,] 0.010498785 0.02099757 0.98950121
[23,] 0.012859546 0.02571909 0.98714045
[24,] 0.032510725 0.06502145 0.96748927
[25,] 0.060702609 0.12140522 0.93929739
[26,] 0.092770865 0.18554173 0.90722913
[27,] 0.145098710 0.29019742 0.85490129
[28,] 0.280143589 0.56028718 0.71985641
[29,] 0.434265911 0.86853182 0.56573409
[30,] 0.533040154 0.93391969 0.46695985
[31,] 0.653948022 0.69210396 0.34605198
[32,] 0.839814832 0.32037034 0.16018517
[33,] 0.921834922 0.15633016 0.07816508
[34,] 0.941601947 0.11679611 0.05839805
[35,] 0.936468972 0.12706206 0.06353103
[36,] 0.938906330 0.12218734 0.06109367
[37,] 0.940684052 0.11863190 0.05931595
[38,] 0.914386337 0.17122733 0.08561366
[39,] 0.878624844 0.24275031 0.12137516
[40,] 0.832583948 0.33483210 0.16741605
[41,] 0.929215627 0.14156875 0.07078437
[42,] 0.919056408 0.16188718 0.08094359
[43,] 0.896799502 0.20640100 0.10320050
[44,] 0.853507478 0.29298504 0.14649252
[45,] 0.887151621 0.22569676 0.11284838
[46,] 0.860999598 0.27800080 0.13900040
[47,] 0.830537519 0.33892496 0.16946248
[48,] 0.864374333 0.27125133 0.13562567
[49,] 0.838851288 0.32229742 0.16114871
[50,] 0.909027649 0.18194470 0.09097235
[51,] 0.915070464 0.16985907 0.08492954
[52,] 0.956935534 0.08612893 0.04306447
[53,] 0.910927243 0.17814551 0.08907276
[54,] 0.830191491 0.33961702 0.16980851
[55,] 0.725405709 0.54918858 0.27459429
> postscript(file="/var/wessaorg/rcomp/tmp/19iyk1324137844.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/2prcd1324137844.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/3hcta1324137844.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/4mjpb1324137844.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/5ve0v1324137844.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 = 66
Frequency = 1
1 2 3 4 5 6
18.2277951 -11.2102503 54.7613998 49.2426247 28.3632748 -16.7768137
7 8 9 10 11 12
33.4856388 10.6284416 -0.3631247 -13.5143778 42.7159555 18.2610383
13 14 15 16 17 18
1.8525799 -37.2726633 16.9396453 -19.9774784 -27.9335628 -55.7514537
19 20 21 22 23 24
10.6398757 -18.4459278 -70.6377474 -15.6731755 -50.5143830 -59.5460855
25 26 27 28 29 30
-82.5397839 -26.5473454 -55.5747798 -67.1374170 -78.9288466 -13.5145808
31 32 33 34 35 36
-18.6920980 -15.6208571 -31.0935141 34.8756764 15.7374345 20.7528118
37 38 39 40 41 42
13.0942372 76.3788890 71.6830360 56.8351523 41.0043547 89.6183385
43 44 45 46 47 48
61.9608787 53.5639862 45.2122428 91.5106181 41.7444808 28.9338530
49 50 51 52 53 54
0.2942366 38.0996348 4.3740570 -10.9708194 -44.9217665 -9.3892761
55 56 57 58 59 60
-47.1720397 -38.6639491 -51.6802511 7.8121844 -2.3437606 9.5651294
61 62 63 64 65 66
-25.1815408 26.2312547 -15.7282539 -22.5933920 -48.6175121 -9.8719284
> postscript(file="/var/wessaorg/rcomp/tmp/6yjn81324137844.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 = 66
Frequency = 1
lag(myerror, k = 1) myerror
0 18.2277951 NA
1 -11.2102503 18.2277951
2 54.7613998 -11.2102503
3 49.2426247 54.7613998
4 28.3632748 49.2426247
5 -16.7768137 28.3632748
6 33.4856388 -16.7768137
7 10.6284416 33.4856388
8 -0.3631247 10.6284416
9 -13.5143778 -0.3631247
10 42.7159555 -13.5143778
11 18.2610383 42.7159555
12 1.8525799 18.2610383
13 -37.2726633 1.8525799
14 16.9396453 -37.2726633
15 -19.9774784 16.9396453
16 -27.9335628 -19.9774784
17 -55.7514537 -27.9335628
18 10.6398757 -55.7514537
19 -18.4459278 10.6398757
20 -70.6377474 -18.4459278
21 -15.6731755 -70.6377474
22 -50.5143830 -15.6731755
23 -59.5460855 -50.5143830
24 -82.5397839 -59.5460855
25 -26.5473454 -82.5397839
26 -55.5747798 -26.5473454
27 -67.1374170 -55.5747798
28 -78.9288466 -67.1374170
29 -13.5145808 -78.9288466
30 -18.6920980 -13.5145808
31 -15.6208571 -18.6920980
32 -31.0935141 -15.6208571
33 34.8756764 -31.0935141
34 15.7374345 34.8756764
35 20.7528118 15.7374345
36 13.0942372 20.7528118
37 76.3788890 13.0942372
38 71.6830360 76.3788890
39 56.8351523 71.6830360
40 41.0043547 56.8351523
41 89.6183385 41.0043547
42 61.9608787 89.6183385
43 53.5639862 61.9608787
44 45.2122428 53.5639862
45 91.5106181 45.2122428
46 41.7444808 91.5106181
47 28.9338530 41.7444808
48 0.2942366 28.9338530
49 38.0996348 0.2942366
50 4.3740570 38.0996348
51 -10.9708194 4.3740570
52 -44.9217665 -10.9708194
53 -9.3892761 -44.9217665
54 -47.1720397 -9.3892761
55 -38.6639491 -47.1720397
56 -51.6802511 -38.6639491
57 7.8121844 -51.6802511
58 -2.3437606 7.8121844
59 9.5651294 -2.3437606
60 -25.1815408 9.5651294
61 26.2312547 -25.1815408
62 -15.7282539 26.2312547
63 -22.5933920 -15.7282539
64 -48.6175121 -22.5933920
65 -9.8719284 -48.6175121
66 NA -9.8719284
> dum1 <- dum[2:length(myerror),]
> dum1
lag(myerror, k = 1) myerror
[1,] -11.2102503 18.2277951
[2,] 54.7613998 -11.2102503
[3,] 49.2426247 54.7613998
[4,] 28.3632748 49.2426247
[5,] -16.7768137 28.3632748
[6,] 33.4856388 -16.7768137
[7,] 10.6284416 33.4856388
[8,] -0.3631247 10.6284416
[9,] -13.5143778 -0.3631247
[10,] 42.7159555 -13.5143778
[11,] 18.2610383 42.7159555
[12,] 1.8525799 18.2610383
[13,] -37.2726633 1.8525799
[14,] 16.9396453 -37.2726633
[15,] -19.9774784 16.9396453
[16,] -27.9335628 -19.9774784
[17,] -55.7514537 -27.9335628
[18,] 10.6398757 -55.7514537
[19,] -18.4459278 10.6398757
[20,] -70.6377474 -18.4459278
[21,] -15.6731755 -70.6377474
[22,] -50.5143830 -15.6731755
[23,] -59.5460855 -50.5143830
[24,] -82.5397839 -59.5460855
[25,] -26.5473454 -82.5397839
[26,] -55.5747798 -26.5473454
[27,] -67.1374170 -55.5747798
[28,] -78.9288466 -67.1374170
[29,] -13.5145808 -78.9288466
[30,] -18.6920980 -13.5145808
[31,] -15.6208571 -18.6920980
[32,] -31.0935141 -15.6208571
[33,] 34.8756764 -31.0935141
[34,] 15.7374345 34.8756764
[35,] 20.7528118 15.7374345
[36,] 13.0942372 20.7528118
[37,] 76.3788890 13.0942372
[38,] 71.6830360 76.3788890
[39,] 56.8351523 71.6830360
[40,] 41.0043547 56.8351523
[41,] 89.6183385 41.0043547
[42,] 61.9608787 89.6183385
[43,] 53.5639862 61.9608787
[44,] 45.2122428 53.5639862
[45,] 91.5106181 45.2122428
[46,] 41.7444808 91.5106181
[47,] 28.9338530 41.7444808
[48,] 0.2942366 28.9338530
[49,] 38.0996348 0.2942366
[50,] 4.3740570 38.0996348
[51,] -10.9708194 4.3740570
[52,] -44.9217665 -10.9708194
[53,] -9.3892761 -44.9217665
[54,] -47.1720397 -9.3892761
[55,] -38.6639491 -47.1720397
[56,] -51.6802511 -38.6639491
[57,] 7.8121844 -51.6802511
[58,] -2.3437606 7.8121844
[59,] 9.5651294 -2.3437606
[60,] -25.1815408 9.5651294
[61,] 26.2312547 -25.1815408
[62,] -15.7282539 26.2312547
[63,] -22.5933920 -15.7282539
[64,] -48.6175121 -22.5933920
[65,] -9.8719284 -48.6175121
> z <- as.data.frame(dum1)
> z
lag(myerror, k = 1) myerror
1 -11.2102503 18.2277951
2 54.7613998 -11.2102503
3 49.2426247 54.7613998
4 28.3632748 49.2426247
5 -16.7768137 28.3632748
6 33.4856388 -16.7768137
7 10.6284416 33.4856388
8 -0.3631247 10.6284416
9 -13.5143778 -0.3631247
10 42.7159555 -13.5143778
11 18.2610383 42.7159555
12 1.8525799 18.2610383
13 -37.2726633 1.8525799
14 16.9396453 -37.2726633
15 -19.9774784 16.9396453
16 -27.9335628 -19.9774784
17 -55.7514537 -27.9335628
18 10.6398757 -55.7514537
19 -18.4459278 10.6398757
20 -70.6377474 -18.4459278
21 -15.6731755 -70.6377474
22 -50.5143830 -15.6731755
23 -59.5460855 -50.5143830
24 -82.5397839 -59.5460855
25 -26.5473454 -82.5397839
26 -55.5747798 -26.5473454
27 -67.1374170 -55.5747798
28 -78.9288466 -67.1374170
29 -13.5145808 -78.9288466
30 -18.6920980 -13.5145808
31 -15.6208571 -18.6920980
32 -31.0935141 -15.6208571
33 34.8756764 -31.0935141
34 15.7374345 34.8756764
35 20.7528118 15.7374345
36 13.0942372 20.7528118
37 76.3788890 13.0942372
38 71.6830360 76.3788890
39 56.8351523 71.6830360
40 41.0043547 56.8351523
41 89.6183385 41.0043547
42 61.9608787 89.6183385
43 53.5639862 61.9608787
44 45.2122428 53.5639862
45 91.5106181 45.2122428
46 41.7444808 91.5106181
47 28.9338530 41.7444808
48 0.2942366 28.9338530
49 38.0996348 0.2942366
50 4.3740570 38.0996348
51 -10.9708194 4.3740570
52 -44.9217665 -10.9708194
53 -9.3892761 -44.9217665
54 -47.1720397 -9.3892761
55 -38.6639491 -47.1720397
56 -51.6802511 -38.6639491
57 7.8121844 -51.6802511
58 -2.3437606 7.8121844
59 9.5651294 -2.3437606
60 -25.1815408 9.5651294
61 26.2312547 -25.1815408
62 -15.7282539 26.2312547
63 -22.5933920 -15.7282539
64 -48.6175121 -22.5933920
65 -9.8719284 -48.6175121
> 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/7jmiv1324137844.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/8938i1324137844.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/9m62b1324137844.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/10xfof1324137844.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/11siqe1324137844.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/122fqy1324137844.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/13uqts1324137844.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/14ah1u1324137844.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/153jpc1324137844.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/16nyv91324137844.tab")
+ }
>
> try(system("convert tmp/19iyk1324137844.ps tmp/19iyk1324137844.png",intern=TRUE))
character(0)
> try(system("convert tmp/2prcd1324137844.ps tmp/2prcd1324137844.png",intern=TRUE))
character(0)
> try(system("convert tmp/3hcta1324137844.ps tmp/3hcta1324137844.png",intern=TRUE))
character(0)
> try(system("convert tmp/4mjpb1324137844.ps tmp/4mjpb1324137844.png",intern=TRUE))
character(0)
> try(system("convert tmp/5ve0v1324137844.ps tmp/5ve0v1324137844.png",intern=TRUE))
character(0)
> try(system("convert tmp/6yjn81324137844.ps tmp/6yjn81324137844.png",intern=TRUE))
character(0)
> try(system("convert tmp/7jmiv1324137844.ps tmp/7jmiv1324137844.png",intern=TRUE))
character(0)
> try(system("convert tmp/8938i1324137844.ps tmp/8938i1324137844.png",intern=TRUE))
character(0)
> try(system("convert tmp/9m62b1324137844.ps tmp/9m62b1324137844.png",intern=TRUE))
character(0)
> try(system("convert tmp/10xfof1324137844.ps tmp/10xfof1324137844.png",intern=TRUE))
character(0)
>
>
> proc.time()
user system elapsed
3.161 0.584 3.779