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(186448
+ ,17822
+ ,1942
+ ,16739
+ ,4872
+ ,1020
+ ,190530
+ ,22422
+ ,2547
+ ,17851
+ ,4905
+ ,1200
+ ,194207
+ ,18817
+ ,2033
+ ,17034
+ ,4971
+ ,1279
+ ,190855
+ ,22043
+ ,2049
+ ,18055
+ ,4971
+ ,1308
+ ,200779
+ ,19191
+ ,2007
+ ,18216
+ ,4930
+ ,1173
+ ,204428
+ ,23171
+ ,2660
+ ,18960
+ ,5001
+ ,1291
+ ,207617
+ ,19463
+ ,2063
+ ,17903
+ ,5059
+ ,1466
+ ,212071
+ ,22522
+ ,2113
+ ,18842
+ ,5085
+ ,1507
+ ,214239
+ ,20265
+ ,2145
+ ,18907
+ ,5111
+ ,1478
+ ,215883
+ ,24249
+ ,2866
+ ,19862
+ ,5190
+ ,1629
+ ,223484
+ ,20299
+ ,2163
+ ,18836
+ ,5076
+ ,1712
+ ,221529
+ ,25455
+ ,2157
+ ,19846
+ ,5134
+ ,1727
+ ,225247
+ ,21089
+ ,2201
+ ,19511
+ ,4804
+ ,1519
+ ,226699
+ ,26237
+ ,2838
+ ,20318
+ ,4579
+ ,1617
+ ,231406
+ ,21362
+ ,2142
+ ,19843
+ ,4526
+ ,1637
+ ,232324
+ ,26489
+ ,2253
+ ,20975
+ ,4550
+ ,1633
+ ,237192
+ ,21828
+ ,2258
+ ,20485
+ ,4566
+ ,1469
+ ,236727
+ ,27496
+ ,2979
+ ,21407
+ ,4588
+ ,1657
+ ,240698
+ ,21991
+ ,2288
+ ,20404
+ ,4564
+ ,1599
+ ,240688
+ ,27611
+ ,2431
+ ,21454
+ ,4723
+ ,1420
+ ,245283
+ ,22512
+ ,2393
+ ,21558
+ ,4553
+ ,1495
+ ,243556
+ ,28581
+ ,3244
+ ,22442
+ ,4556
+ ,1623
+ ,247826
+ ,23000
+ ,2476
+ ,21201
+ ,4542
+ ,1346
+ ,245798
+ ,28385
+ ,2490
+ ,21804
+ ,4234
+ ,1613
+ ,250479
+ ,23387
+ ,2547
+ ,22537
+ ,4341
+ ,1563
+ ,249216
+ ,30192
+ ,3461
+ ,22736
+ ,4269
+ ,2071
+ ,251896
+ ,24346
+ ,2549
+ ,21525
+ ,4217
+ ,1584
+ ,247616
+ ,30393
+ ,2496
+ ,22427
+ ,4207
+ ,1843
+ ,249994
+ ,24753
+ ,2532
+ ,23437
+ ,4267
+ ,1598
+ ,246552
+ ,31723
+ ,3553
+ ,23366
+ ,4249
+ ,1687
+ ,248771
+ ,24838
+ ,2555
+ ,22281
+ ,4217
+ ,1473
+ ,247551
+ ,32272
+ ,2565
+ ,22994
+ ,4172
+ ,2080
+ ,249745
+ ,25219
+ ,2548
+ ,24007
+ ,4161
+ ,1703
+ ,245742
+ ,33191
+ ,3932
+ ,24145
+ ,4103
+ ,1832
+ ,249019
+ ,26218
+ ,2525
+ ,23065
+ ,4027
+ ,1781
+ ,245841
+ ,33537
+ ,2633
+ ,24374
+ ,4042
+ ,2481
+ ,248771
+ ,27975
+ ,2657
+ ,24805
+ ,4120
+ ,1977
+ ,244723
+ ,34356
+ ,3829
+ ,25159
+ ,4188
+ ,1974
+ ,246878
+ ,27082
+ ,2769
+ ,23751
+ ,4185
+ ,1777
+ ,246014
+ ,34333
+ ,2816
+ ,25487
+ ,4216
+ ,2303
+ ,248496
+ ,28141
+ ,3052
+ ,25608
+ ,4250
+ ,1480
+ ,244351
+ ,36125
+ ,4146
+ ,26396
+ ,4259
+ ,1907
+ ,248016
+ ,28451
+ ,3185
+ ,25207
+ ,4206
+ ,1610
+ ,246509
+ ,35801
+ ,3147
+ ,27000
+ ,4132
+ ,1546
+ ,249426
+ ,28979
+ ,3161
+ ,27369
+ ,3944
+ ,1718
+ ,247840
+ ,37285
+ ,4311
+ ,28401
+ ,3872
+ ,1841
+ ,251035
+ ,30310
+ ,3155
+ ,27126
+ ,3797
+ ,1650
+ ,250161
+ ,36721
+ ,3284
+ ,28474
+ ,3840
+ ,1671
+ ,254278
+ ,29534
+ ,3350
+ ,28926
+ ,3895
+ ,1974
+ ,250801
+ ,38626
+ ,4268
+ ,29894
+ ,3633
+ ,2153
+ ,253985
+ ,29654
+ ,3220
+ ,28822
+ ,3622
+ ,1898
+ ,249174
+ ,42638
+ ,8289
+ ,29849
+ ,3562
+ ,2725
+ ,251287
+ ,31372
+ ,3419
+ ,30624
+ ,3555
+ ,2047
+ ,247947
+ ,39603
+ ,3902
+ ,31038
+ ,3489
+ ,1698
+ ,249992
+ ,31647
+ ,3223
+ ,29468
+ ,3500
+ ,1768
+ ,243805
+ ,39946
+ ,3447
+ ,31294
+ ,3373
+ ,1921
+ ,255812
+ ,31518
+ ,3389
+ ,32110
+ ,3285
+ ,9782
+ ,250417
+ ,42743
+ ,4637
+ ,32827
+ ,3292
+ ,2231
+ ,253033
+ ,33462
+ ,3509
+ ,31327
+ ,3241
+ ,2062
+ ,248705
+ ,41744
+ ,4107
+ ,32749
+ ,3266
+ ,2132
+ ,253950
+ ,33142
+ ,3632
+ ,33598
+ ,3168
+ ,2465
+ ,251484
+ ,41753
+ ,4490
+ ,33878
+ ,3181
+ ,2198
+ ,251093
+ ,35487
+ ,3649
+ ,32292
+ ,3246
+ ,2330
+ ,245996
+ ,44720
+ ,3983
+ ,34021
+ ,3159
+ ,1214
+ ,252721
+ ,33472
+ ,3678
+ ,34955
+ ,3209
+ ,2517
+ ,248019
+ ,45134
+ ,4570
+ ,35322
+ ,3220
+ ,2255
+ ,250464
+ ,36255
+ ,3778
+ ,33816
+ ,3305
+ ,2379
+ ,245571
+ ,46228
+ ,4153
+ ,35766
+ ,3251
+ ,2349
+ ,252690
+ ,35483
+ ,4027
+ ,36770
+ ,3281
+ ,2219
+ ,250183
+ ,47663
+ ,5050
+ ,37762
+ ,3304
+ ,2470
+ ,253639
+ ,38064
+ ,4155
+ ,36298
+ ,3270
+ ,2540
+ ,254436
+ ,47177
+ ,4475
+ ,39219
+ ,3377
+ ,2667
+ ,265280
+ ,35062
+ ,4117
+ ,39664
+ ,3235
+ ,3507
+ ,268705
+ ,45062
+ ,5193
+ ,40178
+ ,3125
+ ,2972
+ ,270643
+ ,36943
+ ,4199
+ ,38402
+ ,3091
+ ,2678
+ ,271480
+ ,46194
+ ,4391
+ ,40957
+ ,3102
+ ,2979)
+ ,dim=c(6
+ ,76)
+ ,dimnames=list(c('Nettoschuld'
+ ,'Fiscale_en_parafiscale_ontvangsten'
+ ,'Niet-fiscale_en_niet-parafiscale_ontvangsten'
+ ,'Lopende_uitgaven_exclusief_rentelasten'
+ ,'Rentelasten'
+ ,'Kapitaaluitgaven')
+ ,1:76))
> y <- array(NA,dim=c(6,76),dimnames=list(c('Nettoschuld','Fiscale_en_parafiscale_ontvangsten','Niet-fiscale_en_niet-parafiscale_ontvangsten','Lopende_uitgaven_exclusief_rentelasten','Rentelasten','Kapitaaluitgaven'),1:76))
> 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 = 'Include Quarterly 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
Nettoschuld Fiscale_en_parafiscale_ontvangsten
1 186448 17822
2 190530 22422
3 194207 18817
4 190855 22043
5 200779 19191
6 204428 23171
7 207617 19463
8 212071 22522
9 214239 20265
10 215883 24249
11 223484 20299
12 221529 25455
13 225247 21089
14 226699 26237
15 231406 21362
16 232324 26489
17 237192 21828
18 236727 27496
19 240698 21991
20 240688 27611
21 245283 22512
22 243556 28581
23 247826 23000
24 245798 28385
25 250479 23387
26 249216 30192
27 251896 24346
28 247616 30393
29 249994 24753
30 246552 31723
31 248771 24838
32 247551 32272
33 249745 25219
34 245742 33191
35 249019 26218
36 245841 33537
37 248771 27975
38 244723 34356
39 246878 27082
40 246014 34333
41 248496 28141
42 244351 36125
43 248016 28451
44 246509 35801
45 249426 28979
46 247840 37285
47 251035 30310
48 250161 36721
49 254278 29534
50 250801 38626
51 253985 29654
52 249174 42638
53 251287 31372
54 247947 39603
55 249992 31647
56 243805 39946
57 255812 31518
58 250417 42743
59 253033 33462
60 248705 41744
61 253950 33142
62 251484 41753
63 251093 35487
64 245996 44720
65 252721 33472
66 248019 45134
67 250464 36255
68 245571 46228
69 252690 35483
70 250183 47663
71 253639 38064
72 254436 47177
73 265280 35062
74 268705 45062
75 270643 36943
76 271480 46194
Niet-fiscale_en_niet-parafiscale_ontvangsten
1 1942
2 2547
3 2033
4 2049
5 2007
6 2660
7 2063
8 2113
9 2145
10 2866
11 2163
12 2157
13 2201
14 2838
15 2142
16 2253
17 2258
18 2979
19 2288
20 2431
21 2393
22 3244
23 2476
24 2490
25 2547
26 3461
27 2549
28 2496
29 2532
30 3553
31 2555
32 2565
33 2548
34 3932
35 2525
36 2633
37 2657
38 3829
39 2769
40 2816
41 3052
42 4146
43 3185
44 3147
45 3161
46 4311
47 3155
48 3284
49 3350
50 4268
51 3220
52 8289
53 3419
54 3902
55 3223
56 3447
57 3389
58 4637
59 3509
60 4107
61 3632
62 4490
63 3649
64 3983
65 3678
66 4570
67 3778
68 4153
69 4027
70 5050
71 4155
72 4475
73 4117
74 5193
75 4199
76 4391
Lopende_uitgaven_exclusief_rentelasten Rentelasten Kapitaaluitgaven Q1 Q2 Q3
1 16739 4872 1020 1 0 0
2 17851 4905 1200 0 1 0
3 17034 4971 1279 0 0 1
4 18055 4971 1308 0 0 0
5 18216 4930 1173 1 0 0
6 18960 5001 1291 0 1 0
7 17903 5059 1466 0 0 1
8 18842 5085 1507 0 0 0
9 18907 5111 1478 1 0 0
10 19862 5190 1629 0 1 0
11 18836 5076 1712 0 0 1
12 19846 5134 1727 0 0 0
13 19511 4804 1519 1 0 0
14 20318 4579 1617 0 1 0
15 19843 4526 1637 0 0 1
16 20975 4550 1633 0 0 0
17 20485 4566 1469 1 0 0
18 21407 4588 1657 0 1 0
19 20404 4564 1599 0 0 1
20 21454 4723 1420 0 0 0
21 21558 4553 1495 1 0 0
22 22442 4556 1623 0 1 0
23 21201 4542 1346 0 0 1
24 21804 4234 1613 0 0 0
25 22537 4341 1563 1 0 0
26 22736 4269 2071 0 1 0
27 21525 4217 1584 0 0 1
28 22427 4207 1843 0 0 0
29 23437 4267 1598 1 0 0
30 23366 4249 1687 0 1 0
31 22281 4217 1473 0 0 1
32 22994 4172 2080 0 0 0
33 24007 4161 1703 1 0 0
34 24145 4103 1832 0 1 0
35 23065 4027 1781 0 0 1
36 24374 4042 2481 0 0 0
37 24805 4120 1977 1 0 0
38 25159 4188 1974 0 1 0
39 23751 4185 1777 0 0 1
40 25487 4216 2303 0 0 0
41 25608 4250 1480 1 0 0
42 26396 4259 1907 0 1 0
43 25207 4206 1610 0 0 1
44 27000 4132 1546 0 0 0
45 27369 3944 1718 1 0 0
46 28401 3872 1841 0 1 0
47 27126 3797 1650 0 0 1
48 28474 3840 1671 0 0 0
49 28926 3895 1974 1 0 0
50 29894 3633 2153 0 1 0
51 28822 3622 1898 0 0 1
52 29849 3562 2725 0 0 0
53 30624 3555 2047 1 0 0
54 31038 3489 1698 0 1 0
55 29468 3500 1768 0 0 1
56 31294 3373 1921 0 0 0
57 32110 3285 9782 1 0 0
58 32827 3292 2231 0 1 0
59 31327 3241 2062 0 0 1
60 32749 3266 2132 0 0 0
61 33598 3168 2465 1 0 0
62 33878 3181 2198 0 1 0
63 32292 3246 2330 0 0 1
64 34021 3159 1214 0 0 0
65 34955 3209 2517 1 0 0
66 35322 3220 2255 0 1 0
67 33816 3305 2379 0 0 1
68 35766 3251 2349 0 0 0
69 36770 3281 2219 1 0 0
70 37762 3304 2470 0 1 0
71 36298 3270 2540 0 0 1
72 39219 3377 2667 0 0 0
73 39664 3235 3507 1 0 0
74 40178 3125 2972 0 1 0
75 38402 3091 2678 0 0 1
76 40957 3102 2979 0 0 0
t
1 1
2 2
3 3
4 4
5 5
6 6
7 7
8 8
9 9
10 10
11 11
12 12
13 13
14 14
15 15
16 16
17 17
18 18
19 19
20 20
21 21
22 22
23 23
24 24
25 25
26 26
27 27
28 28
29 29
30 30
31 31
32 32
33 33
34 34
35 35
36 36
37 37
38 38
39 39
40 40
41 41
42 42
43 43
44 44
45 45
46 46
47 47
48 48
49 49
50 50
51 51
52 52
53 53
54 54
55 55
56 56
57 57
58 58
59 59
60 60
61 61
62 62
63 63
64 64
65 65
66 66
67 67
68 68
69 69
70 70
71 71
72 72
73 73
74 74
75 75
76 76
> k <- length(x[1,])
> df <- as.data.frame(x)
> (mylm <- lm(df))
Call:
lm(formula = df)
Coefficients:
(Intercept)
2.983e+05
Fiscale_en_parafiscale_ontvangsten
-1.950e+00
`Niet-fiscale_en_niet-parafiscale_ontvangsten`
7.302e-01
Lopende_uitgaven_exclusief_rentelasten
-4.036e+00
Rentelasten
4.006e+00
Kapitaaluitgaven
7.384e-01
Q1
-1.137e+04
Q2
1.758e+03
Q3
-1.491e+04
t
2.520e+03
> (mysum <- summary(mylm))
Call:
lm(formula = df)
Residuals:
Min 1Q Median 3Q Max
-24041 -6016 -1105 6218 19233
Coefficients:
Estimate Std. Error t value
(Intercept) 2.983e+05 4.583e+04 6.509
Fiscale_en_parafiscale_ontvangsten -1.950e+00 9.609e-01 -2.030
`Niet-fiscale_en_niet-parafiscale_ontvangsten` 7.302e-01 2.270e+00 0.322
Lopende_uitgaven_exclusief_rentelasten -4.035e+00 9.352e-01 -4.315
Rentelasten 4.006e+00 7.532e+00 0.532
Kapitaaluitgaven 7.384e-01 1.324e+00 0.558
Q1 -1.137e+04 7.243e+03 -1.570
Q2 1.758e+03 3.440e+03 0.511
Q3 -1.491e+04 7.091e+03 -2.103
t 2.520e+03 4.371e+02 5.764
Pr(>|t|)
(Intercept) 1.20e-08 ***
Fiscale_en_parafiscale_ontvangsten 0.0464 *
`Niet-fiscale_en_niet-parafiscale_ontvangsten` 0.7488
Lopende_uitgaven_exclusief_rentelasten 5.46e-05 ***
Rentelasten 0.5966
Kapitaaluitgaven 0.5790
Q1 0.1213
Q2 0.6110
Q3 0.0393 *
t 2.36e-07 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 9859 on 66 degrees of freedom
Multiple R-squared: 0.7421, Adjusted R-squared: 0.707
F-statistic: 21.1 on 9 and 66 DF, p-value: 3.075e-16
> if (n > n25) {
+ kp3 <- k + 3
+ nmkm3 <- n - k - 3
+ gqarr <- array(NA, dim=c(nmkm3-kp3+1,3))
+ numgqtests <- 0
+ numsignificant1 <- 0
+ numsignificant5 <- 0
+ numsignificant10 <- 0
+ for (mypoint in kp3:nmkm3) {
+ j <- 0
+ numgqtests <- numgqtests + 1
+ for (myalt in c('greater', 'two.sided', 'less')) {
+ j <- j + 1
+ gqarr[mypoint-kp3+1,j] <- gqtest(mylm, point=mypoint, alternative=myalt)$p.value
+ }
+ if (gqarr[mypoint-kp3+1,2] < 0.01) numsignificant1 <- numsignificant1 + 1
+ if (gqarr[mypoint-kp3+1,2] < 0.05) numsignificant5 <- numsignificant5 + 1
+ if (gqarr[mypoint-kp3+1,2] < 0.10) numsignificant10 <- numsignificant10 + 1
+ }
+ gqarr
+ }
[,1] [,2] [,3]
[1,] 0.13254987 2.650997e-01 8.674501e-01
[2,] 0.08143398 1.628680e-01 9.185660e-01
[3,] 0.04988052 9.976105e-02 9.501195e-01
[4,] 0.06416912 1.283382e-01 9.358309e-01
[5,] 0.03836368 7.672736e-02 9.616363e-01
[6,] 0.03358988 6.717976e-02 9.664101e-01
[7,] 0.03760166 7.520332e-02 9.623983e-01
[8,] 0.04147404 8.294808e-02 9.585260e-01
[9,] 0.03517103 7.034206e-02 9.648290e-01
[10,] 0.04202992 8.405983e-02 9.579701e-01
[11,] 0.04000834 8.001669e-02 9.599917e-01
[12,] 0.06358779 1.271756e-01 9.364122e-01
[13,] 0.10622341 2.124468e-01 8.937766e-01
[14,] 0.07978189 1.595638e-01 9.202181e-01
[15,] 0.15775078 3.155016e-01 8.422492e-01
[16,] 0.70834761 5.833048e-01 2.916524e-01
[17,] 0.99232392 1.535216e-02 7.676078e-03
[18,] 0.99636943 7.261136e-03 3.630568e-03
[19,] 0.99931369 1.372618e-03 6.863090e-04
[20,] 0.99970203 5.959407e-04 2.979704e-04
[21,] 0.99998109 3.781233e-05 1.890616e-05
[22,] 0.99999256 1.488789e-05 7.443946e-06
[23,] 0.99999704 5.913208e-06 2.956604e-06
[24,] 0.99999953 9.405203e-07 4.702601e-07
[25,] 0.99999984 3.131240e-07 1.565620e-07
[26,] 0.99999992 1.648523e-07 8.242616e-08
[27,] 0.99999990 1.942427e-07 9.712137e-08
[28,] 0.99999991 1.843427e-07 9.217134e-08
[29,] 0.99999996 7.936611e-08 3.968306e-08
[30,] 0.99999992 1.601754e-07 8.008770e-08
[31,] 0.99999985 3.023171e-07 1.511586e-07
[32,] 0.99999958 8.493989e-07 4.246994e-07
[33,] 0.99999976 4.743692e-07 2.371846e-07
[34,] 0.99999950 1.006450e-06 5.032250e-07
[35,] 0.99999876 2.484340e-06 1.242170e-06
[36,] 0.99999653 6.948387e-06 3.474193e-06
[37,] 0.99999993 1.409122e-07 7.045608e-08
[38,] 0.99999981 3.715788e-07 1.857894e-07
[39,] 0.99999942 1.165633e-06 5.828164e-07
[40,] 0.99999786 4.275602e-06 2.137801e-06
[41,] 0.99999981 3.837903e-07 1.918951e-07
[42,] 0.99999939 1.217945e-06 6.089725e-07
[43,] 0.99999995 1.082669e-07 5.413344e-08
[44,] 0.99999965 7.020180e-07 3.510090e-07
[45,] 0.99999859 2.819480e-06 1.409740e-06
[46,] 0.99999671 6.573205e-06 3.286603e-06
[47,] 0.99997857 4.285542e-05 2.142771e-05
[48,] 0.99991462 1.707514e-04 8.537568e-05
[49,] 0.99950858 9.828478e-04 4.914239e-04
[50,] 0.99765025 4.699493e-03 2.349747e-03
[51,] 0.99642721 7.145589e-03 3.572795e-03
> postscript(file="/var/www/html/rcomp/tmp/16r471293202145.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/2zila1293202145.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/3zila1293202145.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/4zila1293202145.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/5ar3d1293202145.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 = 76
Frequency = 1
1 2 3 4 5 6
-22375.2343 -21188.5771 -13637.4839 -24041.2535 -9885.3041 -11967.3889
7 8 9 10 11 12
-6051.9383 -9444.1574 -2671.8803 -6005.7355 4809.1713 -683.8419
13 14 15 16 17 18
3460.5222 2925.7936 11063.2948 8943.8538 11647.3548 9556.6443
19 20 21 22 23 24
13536.4822 10684.7030 15258.0603 12559.5655 15907.6059 10411.1952
25 26 27 28 29 30
16719.2020 13129.6390 14904.3457 8514.7538 12733.1012 6211.7057
31 32 33 34 35 36
5788.6140 4238.6122 5949.8710 1530.7470 2368.7349 661.0381
37 38 39 40 41 42
3374.8543 -3572.9940 -6205.6751 -3899.8968 -3856.7637 -6048.0905
43 44 45 46 47 48
-6865.2692 -3861.0488 -3294.1405 -805.7459 -924.4814 -1567.7778
49 50 51 52 53 54
-1285.8913 1475.3145 -2017.9334 1136.1443 -2660.9618 -3755.2723
55 56 57 58 59 60
-9013.2415 -8844.2442 -6541.3382 1838.3944 -4397.5556 -4853.7331
61 62 63 64 65 66
-3615.6735 -4286.4464 -8892.7831 -5507.1479 -9039.8166 -5665.6358
67 68 69 70 71 72
-12319.2558 -7358.5972 -8226.2102 352.8943 -5932.6805 6238.4963
73 74 75 76
4310.2485 13715.1882 7880.0490 19232.9020
> postscript(file="/var/www/html/rcomp/tmp/6ar3d1293202145.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 = 76
Frequency = 1
lag(myerror, k = 1) myerror
0 -22375.2343 NA
1 -21188.5771 -22375.2343
2 -13637.4839 -21188.5771
3 -24041.2535 -13637.4839
4 -9885.3041 -24041.2535
5 -11967.3889 -9885.3041
6 -6051.9383 -11967.3889
7 -9444.1574 -6051.9383
8 -2671.8803 -9444.1574
9 -6005.7355 -2671.8803
10 4809.1713 -6005.7355
11 -683.8419 4809.1713
12 3460.5222 -683.8419
13 2925.7936 3460.5222
14 11063.2948 2925.7936
15 8943.8538 11063.2948
16 11647.3548 8943.8538
17 9556.6443 11647.3548
18 13536.4822 9556.6443
19 10684.7030 13536.4822
20 15258.0603 10684.7030
21 12559.5655 15258.0603
22 15907.6059 12559.5655
23 10411.1952 15907.6059
24 16719.2020 10411.1952
25 13129.6390 16719.2020
26 14904.3457 13129.6390
27 8514.7538 14904.3457
28 12733.1012 8514.7538
29 6211.7057 12733.1012
30 5788.6140 6211.7057
31 4238.6122 5788.6140
32 5949.8710 4238.6122
33 1530.7470 5949.8710
34 2368.7349 1530.7470
35 661.0381 2368.7349
36 3374.8543 661.0381
37 -3572.9940 3374.8543
38 -6205.6751 -3572.9940
39 -3899.8968 -6205.6751
40 -3856.7637 -3899.8968
41 -6048.0905 -3856.7637
42 -6865.2692 -6048.0905
43 -3861.0488 -6865.2692
44 -3294.1405 -3861.0488
45 -805.7459 -3294.1405
46 -924.4814 -805.7459
47 -1567.7778 -924.4814
48 -1285.8913 -1567.7778
49 1475.3145 -1285.8913
50 -2017.9334 1475.3145
51 1136.1443 -2017.9334
52 -2660.9618 1136.1443
53 -3755.2723 -2660.9618
54 -9013.2415 -3755.2723
55 -8844.2442 -9013.2415
56 -6541.3382 -8844.2442
57 1838.3944 -6541.3382
58 -4397.5556 1838.3944
59 -4853.7331 -4397.5556
60 -3615.6735 -4853.7331
61 -4286.4464 -3615.6735
62 -8892.7831 -4286.4464
63 -5507.1479 -8892.7831
64 -9039.8166 -5507.1479
65 -5665.6358 -9039.8166
66 -12319.2558 -5665.6358
67 -7358.5972 -12319.2558
68 -8226.2102 -7358.5972
69 352.8943 -8226.2102
70 -5932.6805 352.8943
71 6238.4963 -5932.6805
72 4310.2485 6238.4963
73 13715.1882 4310.2485
74 7880.0490 13715.1882
75 19232.9020 7880.0490
76 NA 19232.9020
> dum1 <- dum[2:length(myerror),]
> dum1
lag(myerror, k = 1) myerror
[1,] -21188.5771 -22375.2343
[2,] -13637.4839 -21188.5771
[3,] -24041.2535 -13637.4839
[4,] -9885.3041 -24041.2535
[5,] -11967.3889 -9885.3041
[6,] -6051.9383 -11967.3889
[7,] -9444.1574 -6051.9383
[8,] -2671.8803 -9444.1574
[9,] -6005.7355 -2671.8803
[10,] 4809.1713 -6005.7355
[11,] -683.8419 4809.1713
[12,] 3460.5222 -683.8419
[13,] 2925.7936 3460.5222
[14,] 11063.2948 2925.7936
[15,] 8943.8538 11063.2948
[16,] 11647.3548 8943.8538
[17,] 9556.6443 11647.3548
[18,] 13536.4822 9556.6443
[19,] 10684.7030 13536.4822
[20,] 15258.0603 10684.7030
[21,] 12559.5655 15258.0603
[22,] 15907.6059 12559.5655
[23,] 10411.1952 15907.6059
[24,] 16719.2020 10411.1952
[25,] 13129.6390 16719.2020
[26,] 14904.3457 13129.6390
[27,] 8514.7538 14904.3457
[28,] 12733.1012 8514.7538
[29,] 6211.7057 12733.1012
[30,] 5788.6140 6211.7057
[31,] 4238.6122 5788.6140
[32,] 5949.8710 4238.6122
[33,] 1530.7470 5949.8710
[34,] 2368.7349 1530.7470
[35,] 661.0381 2368.7349
[36,] 3374.8543 661.0381
[37,] -3572.9940 3374.8543
[38,] -6205.6751 -3572.9940
[39,] -3899.8968 -6205.6751
[40,] -3856.7637 -3899.8968
[41,] -6048.0905 -3856.7637
[42,] -6865.2692 -6048.0905
[43,] -3861.0488 -6865.2692
[44,] -3294.1405 -3861.0488
[45,] -805.7459 -3294.1405
[46,] -924.4814 -805.7459
[47,] -1567.7778 -924.4814
[48,] -1285.8913 -1567.7778
[49,] 1475.3145 -1285.8913
[50,] -2017.9334 1475.3145
[51,] 1136.1443 -2017.9334
[52,] -2660.9618 1136.1443
[53,] -3755.2723 -2660.9618
[54,] -9013.2415 -3755.2723
[55,] -8844.2442 -9013.2415
[56,] -6541.3382 -8844.2442
[57,] 1838.3944 -6541.3382
[58,] -4397.5556 1838.3944
[59,] -4853.7331 -4397.5556
[60,] -3615.6735 -4853.7331
[61,] -4286.4464 -3615.6735
[62,] -8892.7831 -4286.4464
[63,] -5507.1479 -8892.7831
[64,] -9039.8166 -5507.1479
[65,] -5665.6358 -9039.8166
[66,] -12319.2558 -5665.6358
[67,] -7358.5972 -12319.2558
[68,] -8226.2102 -7358.5972
[69,] 352.8943 -8226.2102
[70,] -5932.6805 352.8943
[71,] 6238.4963 -5932.6805
[72,] 4310.2485 6238.4963
[73,] 13715.1882 4310.2485
[74,] 7880.0490 13715.1882
[75,] 19232.9020 7880.0490
> z <- as.data.frame(dum1)
> z
lag(myerror, k = 1) myerror
1 -21188.5771 -22375.2343
2 -13637.4839 -21188.5771
3 -24041.2535 -13637.4839
4 -9885.3041 -24041.2535
5 -11967.3889 -9885.3041
6 -6051.9383 -11967.3889
7 -9444.1574 -6051.9383
8 -2671.8803 -9444.1574
9 -6005.7355 -2671.8803
10 4809.1713 -6005.7355
11 -683.8419 4809.1713
12 3460.5222 -683.8419
13 2925.7936 3460.5222
14 11063.2948 2925.7936
15 8943.8538 11063.2948
16 11647.3548 8943.8538
17 9556.6443 11647.3548
18 13536.4822 9556.6443
19 10684.7030 13536.4822
20 15258.0603 10684.7030
21 12559.5655 15258.0603
22 15907.6059 12559.5655
23 10411.1952 15907.6059
24 16719.2020 10411.1952
25 13129.6390 16719.2020
26 14904.3457 13129.6390
27 8514.7538 14904.3457
28 12733.1012 8514.7538
29 6211.7057 12733.1012
30 5788.6140 6211.7057
31 4238.6122 5788.6140
32 5949.8710 4238.6122
33 1530.7470 5949.8710
34 2368.7349 1530.7470
35 661.0381 2368.7349
36 3374.8543 661.0381
37 -3572.9940 3374.8543
38 -6205.6751 -3572.9940
39 -3899.8968 -6205.6751
40 -3856.7637 -3899.8968
41 -6048.0905 -3856.7637
42 -6865.2692 -6048.0905
43 -3861.0488 -6865.2692
44 -3294.1405 -3861.0488
45 -805.7459 -3294.1405
46 -924.4814 -805.7459
47 -1567.7778 -924.4814
48 -1285.8913 -1567.7778
49 1475.3145 -1285.8913
50 -2017.9334 1475.3145
51 1136.1443 -2017.9334
52 -2660.9618 1136.1443
53 -3755.2723 -2660.9618
54 -9013.2415 -3755.2723
55 -8844.2442 -9013.2415
56 -6541.3382 -8844.2442
57 1838.3944 -6541.3382
58 -4397.5556 1838.3944
59 -4853.7331 -4397.5556
60 -3615.6735 -4853.7331
61 -4286.4464 -3615.6735
62 -8892.7831 -4286.4464
63 -5507.1479 -8892.7831
64 -9039.8166 -5507.1479
65 -5665.6358 -9039.8166
66 -12319.2558 -5665.6358
67 -7358.5972 -12319.2558
68 -8226.2102 -7358.5972
69 352.8943 -8226.2102
70 -5932.6805 352.8943
71 6238.4963 -5932.6805
72 4310.2485 6238.4963
73 13715.1882 4310.2485
74 7880.0490 13715.1882
75 19232.9020 7880.0490
> 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/731kg1293202145.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/831kg1293202145.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/9vs111293202145.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/10vs111293202145.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/11hai71293202145.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/12kbgc1293202145.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/13glel1293202145.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/14k3dr1293202145.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/15n4bx1293202145.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/16qmrl1293202145.tab")
+ }
>
> try(system("convert tmp/16r471293202145.ps tmp/16r471293202145.png",intern=TRUE))
character(0)
> try(system("convert tmp/2zila1293202145.ps tmp/2zila1293202145.png",intern=TRUE))
character(0)
> try(system("convert tmp/3zila1293202145.ps tmp/3zila1293202145.png",intern=TRUE))
character(0)
> try(system("convert tmp/4zila1293202145.ps tmp/4zila1293202145.png",intern=TRUE))
character(0)
> try(system("convert tmp/5ar3d1293202145.ps tmp/5ar3d1293202145.png",intern=TRUE))
character(0)
> try(system("convert tmp/6ar3d1293202145.ps tmp/6ar3d1293202145.png",intern=TRUE))
character(0)
> try(system("convert tmp/731kg1293202145.ps tmp/731kg1293202145.png",intern=TRUE))
character(0)
> try(system("convert tmp/831kg1293202145.ps tmp/831kg1293202145.png",intern=TRUE))
character(0)
> try(system("convert tmp/9vs111293202145.ps tmp/9vs111293202145.png",intern=TRUE))
character(0)
> try(system("convert tmp/10vs111293202145.ps tmp/10vs111293202145.png",intern=TRUE))
character(0)
>
>
> proc.time()
user system elapsed
2.729 1.680 10.931