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(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'
+ ,'ParafiscaleOntvangsten'
+ ,'Niet-parafiscaleOntvangsten'
+ ,'LopendeUitgaven'
+ ,'Rentelasten'
+ ,'KapitaalUitgaven')
+ ,1:76))
> y <- array(NA,dim=c(6,76),dimnames=list(c('Nettoschuld','ParafiscaleOntvangsten','Niet-parafiscaleOntvangsten','LopendeUitgaven','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 = '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
Nettoschuld ParafiscaleOntvangsten Niet-parafiscaleOntvangsten
1 186448 17822 1942
2 190530 22422 2547
3 194207 18817 2033
4 190855 22043 2049
5 200779 19191 2007
6 204428 23171 2660
7 207617 19463 2063
8 212071 22522 2113
9 214239 20265 2145
10 215883 24249 2866
11 223484 20299 2163
12 221529 25455 2157
13 225247 21089 2201
14 226699 26237 2838
15 231406 21362 2142
16 232324 26489 2253
17 237192 21828 2258
18 236727 27496 2979
19 240698 21991 2288
20 240688 27611 2431
21 245283 22512 2393
22 243556 28581 3244
23 247826 23000 2476
24 245798 28385 2490
25 250479 23387 2547
26 249216 30192 3461
27 251896 24346 2549
28 247616 30393 2496
29 249994 24753 2532
30 246552 31723 3553
31 248771 24838 2555
32 247551 32272 2565
33 249745 25219 2548
34 245742 33191 3932
35 249019 26218 2525
36 245841 33537 2633
37 248771 27975 2657
38 244723 34356 3829
39 246878 27082 2769
40 246014 34333 2816
41 248496 28141 3052
42 244351 36125 4146
43 248016 28451 3185
44 246509 35801 3147
45 249426 28979 3161
46 247840 37285 4311
47 251035 30310 3155
48 250161 36721 3284
49 254278 29534 3350
50 250801 38626 4268
51 253985 29654 3220
52 249174 42638 8289
53 251287 31372 3419
54 247947 39603 3902
55 249992 31647 3223
56 243805 39946 3447
57 255812 31518 3389
58 250417 42743 4637
59 253033 33462 3509
60 248705 41744 4107
61 253950 33142 3632
62 251484 41753 4490
63 251093 35487 3649
64 245996 44720 3983
65 252721 33472 3678
66 248019 45134 4570
67 250464 36255 3778
68 245571 46228 4153
69 252690 35483 4027
70 250183 47663 5050
71 253639 38064 4155
72 254436 47177 4475
73 265280 35062 4117
74 268705 45062 5193
75 270643 36943 4199
76 271480 46194 4391
LopendeUitgaven Rentelasten KapitaalUitgaven
1 16739 4872 1020
2 17851 4905 1200
3 17034 4971 1279
4 18055 4971 1308
5 18216 4930 1173
6 18960 5001 1291
7 17903 5059 1466
8 18842 5085 1507
9 18907 5111 1478
10 19862 5190 1629
11 18836 5076 1712
12 19846 5134 1727
13 19511 4804 1519
14 20318 4579 1617
15 19843 4526 1637
16 20975 4550 1633
17 20485 4566 1469
18 21407 4588 1657
19 20404 4564 1599
20 21454 4723 1420
21 21558 4553 1495
22 22442 4556 1623
23 21201 4542 1346
24 21804 4234 1613
25 22537 4341 1563
26 22736 4269 2071
27 21525 4217 1584
28 22427 4207 1843
29 23437 4267 1598
30 23366 4249 1687
31 22281 4217 1473
32 22994 4172 2080
33 24007 4161 1703
34 24145 4103 1832
35 23065 4027 1781
36 24374 4042 2481
37 24805 4120 1977
38 25159 4188 1974
39 23751 4185 1777
40 25487 4216 2303
41 25608 4250 1480
42 26396 4259 1907
43 25207 4206 1610
44 27000 4132 1546
45 27369 3944 1718
46 28401 3872 1841
47 27126 3797 1650
48 28474 3840 1671
49 28926 3895 1974
50 29894 3633 2153
51 28822 3622 1898
52 29849 3562 2725
53 30624 3555 2047
54 31038 3489 1698
55 29468 3500 1768
56 31294 3373 1921
57 32110 3285 9782
58 32827 3292 2231
59 31327 3241 2062
60 32749 3266 2132
61 33598 3168 2465
62 33878 3181 2198
63 32292 3246 2330
64 34021 3159 1214
65 34955 3209 2517
66 35322 3220 2255
67 33816 3305 2379
68 35766 3251 2349
69 36770 3281 2219
70 37762 3304 2470
71 36298 3270 2540
72 39219 3377 2667
73 39664 3235 3507
74 40178 3125 2972
75 38402 3091 2678
76 40957 3102 2979
> k <- length(x[1,])
> df <- as.data.frame(x)
> (mylm <- lm(df))
Call:
lm(formula = df)
Coefficients:
(Intercept) ParafiscaleOntvangsten
3.396e+05 -4.816e-02
`Niet-parafiscaleOntvangsten` LopendeUitgaven
-4.176e-01 -7.901e-02
Rentelasten KapitaalUitgaven
-2.326e+01 5.999e-01
> (mysum <- summary(mylm))
Call:
lm(formula = df)
Residuals:
Min 1Q Median 3Q Max
-37402 -6944 3270 8243 16929
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 3.396e+05 4.329e+04 7.843 3.55e-11 ***
ParafiscaleOntvangsten -4.816e-02 4.515e-01 -0.107 0.915350
`Niet-parafiscaleOntvangsten` -4.176e-01 2.486e+00 -0.168 0.867079
LopendeUitgaven -7.901e-02 7.286e-01 -0.108 0.913949
Rentelasten -2.326e+01 6.586e+00 -3.532 0.000736 ***
KapitaalUitgaven 5.999e-01 1.551e+00 0.387 0.700122
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 11870 on 70 degrees of freedom
Multiple R-squared: 0.6033, Adjusted R-squared: 0.575
F-statistic: 21.29 on 5 and 70 DF, p-value: 6.99e-13
> 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.15260220 3.052044e-01 8.473978e-01
[2,] 0.16160512 3.232102e-01 8.383949e-01
[3,] 0.09254436 1.850887e-01 9.074556e-01
[4,] 0.08599025 1.719805e-01 9.140097e-01
[5,] 0.06852425 1.370485e-01 9.314758e-01
[6,] 0.05728165 1.145633e-01 9.427184e-01
[7,] 0.08596368 1.719274e-01 9.140363e-01
[8,] 0.07380452 1.476090e-01 9.261955e-01
[9,] 0.09447249 1.889450e-01 9.055275e-01
[10,] 0.08225578 1.645116e-01 9.177442e-01
[11,] 0.07913443 1.582689e-01 9.208656e-01
[12,] 0.47406875 9.481375e-01 5.259313e-01
[13,] 0.50162516 9.967497e-01 4.983748e-01
[14,] 0.46525590 9.305118e-01 5.347441e-01
[15,] 0.62330773 7.533845e-01 3.766923e-01
[16,] 0.55436611 8.912678e-01 4.456339e-01
[17,] 0.82952318 3.409536e-01 1.704768e-01
[18,] 0.83104828 3.379034e-01 1.689517e-01
[19,] 0.85729440 2.854112e-01 1.427056e-01
[20,] 0.85805549 2.838890e-01 1.419445e-01
[21,] 0.98963566 2.072867e-02 1.036434e-02
[22,] 0.98506755 2.986490e-02 1.493245e-02
[23,] 0.97829199 4.341602e-02 2.170801e-02
[24,] 0.98759936 2.480128e-02 1.240064e-02
[25,] 0.99917016 1.659677e-03 8.298385e-04
[26,] 0.99915969 1.680622e-03 8.403110e-04
[27,] 0.99935634 1.287320e-03 6.436599e-04
[28,] 0.99991458 1.708493e-04 8.542463e-05
[29,] 0.99997892 4.216548e-05 2.108274e-05
[30,] 0.99998100 3.799200e-05 1.899600e-05
[31,] 0.99997001 5.998685e-05 2.999342e-05
[32,] 0.99997481 5.037882e-05 2.518941e-05
[33,] 0.99997108 5.783157e-05 2.891578e-05
[34,] 0.99996110 7.780459e-05 3.890230e-05
[35,] 0.99992986 1.402773e-04 7.013864e-05
[36,] 0.99986756 2.648724e-04 1.324362e-04
[37,] 0.99990097 1.980518e-04 9.902589e-05
[38,] 0.99986703 2.659330e-04 1.329665e-04
[39,] 0.99981947 3.610684e-04 1.805342e-04
[40,] 0.99975388 4.922497e-04 2.461249e-04
[41,] 0.99972931 5.413733e-04 2.706867e-04
[42,] 0.99980452 3.909516e-04 1.954758e-04
[43,] 0.99989104 2.179211e-04 1.089606e-04
[44,] 0.99975940 4.811908e-04 2.405954e-04
[45,] 0.99975916 4.816893e-04 2.408447e-04
[46,] 0.99977790 4.442029e-04 2.221014e-04
[47,] 0.99996714 6.571411e-05 3.285706e-05
[48,] 0.99997458 5.083677e-05 2.541839e-05
[49,] 0.99998567 2.865800e-05 1.432900e-05
[50,] 0.99998026 3.948002e-05 1.974001e-05
[51,] 0.99998677 2.645640e-05 1.322820e-05
[52,] 0.99998183 3.633269e-05 1.816635e-05
[53,] 0.99993376 1.324843e-04 6.624216e-05
[54,] 0.99975836 4.832874e-04 2.416437e-04
[55,] 0.99972546 5.490787e-04 2.745393e-04
[56,] 0.99892181 2.156377e-03 1.078188e-03
[57,] 0.99648571 7.028586e-03 3.514293e-03
[58,] 0.98732333 2.535334e-02 1.267667e-02
[59,] 0.98884718 2.230564e-02 1.115282e-02
> postscript(file="/var/www/rcomp/tmp/1xyen1293031357.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/2q8e81293031357.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/3q8e81293031357.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/4q8e81293031357.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/5jzdb1293031357.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
-3.740236e+04 -3.209869e+04 -2.738669e+04 -3.051336e+04 -2.160424e+04
6 7 8 9 10
-1.585136e+04 -1.192962e+04 -6.653038e+03 -3.953063e+03 6.365707e+00
11 12 13 14 15
4.340975e+03 4.051722e+03 8.055137e-02 -3.262680e+03 -3.634590e+02
16 17 18 19 20
1.497921e+03 6.575368e+03 7.156233e+03 9.970835e+03 1.418003e+04
21 22 23 24 25
1.452248e+04 1.350598e+04 1.692895e+04 7.889303e+03 1.493021e+04
26 27 28 29 30
1.241281e+04 1.341733e+04 9.089722e+03 1.283355e+04 9.675899e+03
31 32 33 34 35
1.044486e+04 8.232526e+03 1.013008e+04 5.673358e+03 6.204425e+03
36 37 38 39 40
3.456422e+03 8.279320e+03 6.639554e+03 7.938732e+03 7.986278e+03
41 42 43 44 45
1.156277e+04 8.274577e+03 1.002009e+04 7.309981e+03 5.457214e+03
46 47 48 49 50
3.084448e+03 3.730073e+03 4.312833e+03 9.244538e+03 4.635490e+02
51 52 53 54 55
2.590218e+03 -1.289335e+03 -1.447415e+03 -5.482430e+03 -4.014324e+03
56 57 58 59 60
-1.260971e+04 -7.731240e+03 -7.315045e+03 -6.820506e+03 -9.848030e+03
61 62 63 64 65
-7.627914e+03 -8.836213e+03 -8.572741e+03 -1.430316e+04 -7.792093e+03
66 67 68 69 70
-1.111790e+04 -7.647476e+03 -1.298757e+04 -5.583541e+03 -6.613939e+03
71 72 73 74 75
-4.942520e+03 -9.294798e+02 5.409749e+03 7.568573e+03 7.945646e+03
76
9.585543e+03
> postscript(file="/var/www/rcomp/tmp/6jzdb1293031357.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 -3.740236e+04 NA
1 -3.209869e+04 -3.740236e+04
2 -2.738669e+04 -3.209869e+04
3 -3.051336e+04 -2.738669e+04
4 -2.160424e+04 -3.051336e+04
5 -1.585136e+04 -2.160424e+04
6 -1.192962e+04 -1.585136e+04
7 -6.653038e+03 -1.192962e+04
8 -3.953063e+03 -6.653038e+03
9 6.365707e+00 -3.953063e+03
10 4.340975e+03 6.365707e+00
11 4.051722e+03 4.340975e+03
12 8.055137e-02 4.051722e+03
13 -3.262680e+03 8.055137e-02
14 -3.634590e+02 -3.262680e+03
15 1.497921e+03 -3.634590e+02
16 6.575368e+03 1.497921e+03
17 7.156233e+03 6.575368e+03
18 9.970835e+03 7.156233e+03
19 1.418003e+04 9.970835e+03
20 1.452248e+04 1.418003e+04
21 1.350598e+04 1.452248e+04
22 1.692895e+04 1.350598e+04
23 7.889303e+03 1.692895e+04
24 1.493021e+04 7.889303e+03
25 1.241281e+04 1.493021e+04
26 1.341733e+04 1.241281e+04
27 9.089722e+03 1.341733e+04
28 1.283355e+04 9.089722e+03
29 9.675899e+03 1.283355e+04
30 1.044486e+04 9.675899e+03
31 8.232526e+03 1.044486e+04
32 1.013008e+04 8.232526e+03
33 5.673358e+03 1.013008e+04
34 6.204425e+03 5.673358e+03
35 3.456422e+03 6.204425e+03
36 8.279320e+03 3.456422e+03
37 6.639554e+03 8.279320e+03
38 7.938732e+03 6.639554e+03
39 7.986278e+03 7.938732e+03
40 1.156277e+04 7.986278e+03
41 8.274577e+03 1.156277e+04
42 1.002009e+04 8.274577e+03
43 7.309981e+03 1.002009e+04
44 5.457214e+03 7.309981e+03
45 3.084448e+03 5.457214e+03
46 3.730073e+03 3.084448e+03
47 4.312833e+03 3.730073e+03
48 9.244538e+03 4.312833e+03
49 4.635490e+02 9.244538e+03
50 2.590218e+03 4.635490e+02
51 -1.289335e+03 2.590218e+03
52 -1.447415e+03 -1.289335e+03
53 -5.482430e+03 -1.447415e+03
54 -4.014324e+03 -5.482430e+03
55 -1.260971e+04 -4.014324e+03
56 -7.731240e+03 -1.260971e+04
57 -7.315045e+03 -7.731240e+03
58 -6.820506e+03 -7.315045e+03
59 -9.848030e+03 -6.820506e+03
60 -7.627914e+03 -9.848030e+03
61 -8.836213e+03 -7.627914e+03
62 -8.572741e+03 -8.836213e+03
63 -1.430316e+04 -8.572741e+03
64 -7.792093e+03 -1.430316e+04
65 -1.111790e+04 -7.792093e+03
66 -7.647476e+03 -1.111790e+04
67 -1.298757e+04 -7.647476e+03
68 -5.583541e+03 -1.298757e+04
69 -6.613939e+03 -5.583541e+03
70 -4.942520e+03 -6.613939e+03
71 -9.294798e+02 -4.942520e+03
72 5.409749e+03 -9.294798e+02
73 7.568573e+03 5.409749e+03
74 7.945646e+03 7.568573e+03
75 9.585543e+03 7.945646e+03
76 NA 9.585543e+03
> dum1 <- dum[2:length(myerror),]
> dum1
lag(myerror, k = 1) myerror
[1,] -3.209869e+04 -3.740236e+04
[2,] -2.738669e+04 -3.209869e+04
[3,] -3.051336e+04 -2.738669e+04
[4,] -2.160424e+04 -3.051336e+04
[5,] -1.585136e+04 -2.160424e+04
[6,] -1.192962e+04 -1.585136e+04
[7,] -6.653038e+03 -1.192962e+04
[8,] -3.953063e+03 -6.653038e+03
[9,] 6.365707e+00 -3.953063e+03
[10,] 4.340975e+03 6.365707e+00
[11,] 4.051722e+03 4.340975e+03
[12,] 8.055137e-02 4.051722e+03
[13,] -3.262680e+03 8.055137e-02
[14,] -3.634590e+02 -3.262680e+03
[15,] 1.497921e+03 -3.634590e+02
[16,] 6.575368e+03 1.497921e+03
[17,] 7.156233e+03 6.575368e+03
[18,] 9.970835e+03 7.156233e+03
[19,] 1.418003e+04 9.970835e+03
[20,] 1.452248e+04 1.418003e+04
[21,] 1.350598e+04 1.452248e+04
[22,] 1.692895e+04 1.350598e+04
[23,] 7.889303e+03 1.692895e+04
[24,] 1.493021e+04 7.889303e+03
[25,] 1.241281e+04 1.493021e+04
[26,] 1.341733e+04 1.241281e+04
[27,] 9.089722e+03 1.341733e+04
[28,] 1.283355e+04 9.089722e+03
[29,] 9.675899e+03 1.283355e+04
[30,] 1.044486e+04 9.675899e+03
[31,] 8.232526e+03 1.044486e+04
[32,] 1.013008e+04 8.232526e+03
[33,] 5.673358e+03 1.013008e+04
[34,] 6.204425e+03 5.673358e+03
[35,] 3.456422e+03 6.204425e+03
[36,] 8.279320e+03 3.456422e+03
[37,] 6.639554e+03 8.279320e+03
[38,] 7.938732e+03 6.639554e+03
[39,] 7.986278e+03 7.938732e+03
[40,] 1.156277e+04 7.986278e+03
[41,] 8.274577e+03 1.156277e+04
[42,] 1.002009e+04 8.274577e+03
[43,] 7.309981e+03 1.002009e+04
[44,] 5.457214e+03 7.309981e+03
[45,] 3.084448e+03 5.457214e+03
[46,] 3.730073e+03 3.084448e+03
[47,] 4.312833e+03 3.730073e+03
[48,] 9.244538e+03 4.312833e+03
[49,] 4.635490e+02 9.244538e+03
[50,] 2.590218e+03 4.635490e+02
[51,] -1.289335e+03 2.590218e+03
[52,] -1.447415e+03 -1.289335e+03
[53,] -5.482430e+03 -1.447415e+03
[54,] -4.014324e+03 -5.482430e+03
[55,] -1.260971e+04 -4.014324e+03
[56,] -7.731240e+03 -1.260971e+04
[57,] -7.315045e+03 -7.731240e+03
[58,] -6.820506e+03 -7.315045e+03
[59,] -9.848030e+03 -6.820506e+03
[60,] -7.627914e+03 -9.848030e+03
[61,] -8.836213e+03 -7.627914e+03
[62,] -8.572741e+03 -8.836213e+03
[63,] -1.430316e+04 -8.572741e+03
[64,] -7.792093e+03 -1.430316e+04
[65,] -1.111790e+04 -7.792093e+03
[66,] -7.647476e+03 -1.111790e+04
[67,] -1.298757e+04 -7.647476e+03
[68,] -5.583541e+03 -1.298757e+04
[69,] -6.613939e+03 -5.583541e+03
[70,] -4.942520e+03 -6.613939e+03
[71,] -9.294798e+02 -4.942520e+03
[72,] 5.409749e+03 -9.294798e+02
[73,] 7.568573e+03 5.409749e+03
[74,] 7.945646e+03 7.568573e+03
[75,] 9.585543e+03 7.945646e+03
> z <- as.data.frame(dum1)
> z
lag(myerror, k = 1) myerror
1 -3.209869e+04 -3.740236e+04
2 -2.738669e+04 -3.209869e+04
3 -3.051336e+04 -2.738669e+04
4 -2.160424e+04 -3.051336e+04
5 -1.585136e+04 -2.160424e+04
6 -1.192962e+04 -1.585136e+04
7 -6.653038e+03 -1.192962e+04
8 -3.953063e+03 -6.653038e+03
9 6.365707e+00 -3.953063e+03
10 4.340975e+03 6.365707e+00
11 4.051722e+03 4.340975e+03
12 8.055137e-02 4.051722e+03
13 -3.262680e+03 8.055137e-02
14 -3.634590e+02 -3.262680e+03
15 1.497921e+03 -3.634590e+02
16 6.575368e+03 1.497921e+03
17 7.156233e+03 6.575368e+03
18 9.970835e+03 7.156233e+03
19 1.418003e+04 9.970835e+03
20 1.452248e+04 1.418003e+04
21 1.350598e+04 1.452248e+04
22 1.692895e+04 1.350598e+04
23 7.889303e+03 1.692895e+04
24 1.493021e+04 7.889303e+03
25 1.241281e+04 1.493021e+04
26 1.341733e+04 1.241281e+04
27 9.089722e+03 1.341733e+04
28 1.283355e+04 9.089722e+03
29 9.675899e+03 1.283355e+04
30 1.044486e+04 9.675899e+03
31 8.232526e+03 1.044486e+04
32 1.013008e+04 8.232526e+03
33 5.673358e+03 1.013008e+04
34 6.204425e+03 5.673358e+03
35 3.456422e+03 6.204425e+03
36 8.279320e+03 3.456422e+03
37 6.639554e+03 8.279320e+03
38 7.938732e+03 6.639554e+03
39 7.986278e+03 7.938732e+03
40 1.156277e+04 7.986278e+03
41 8.274577e+03 1.156277e+04
42 1.002009e+04 8.274577e+03
43 7.309981e+03 1.002009e+04
44 5.457214e+03 7.309981e+03
45 3.084448e+03 5.457214e+03
46 3.730073e+03 3.084448e+03
47 4.312833e+03 3.730073e+03
48 9.244538e+03 4.312833e+03
49 4.635490e+02 9.244538e+03
50 2.590218e+03 4.635490e+02
51 -1.289335e+03 2.590218e+03
52 -1.447415e+03 -1.289335e+03
53 -5.482430e+03 -1.447415e+03
54 -4.014324e+03 -5.482430e+03
55 -1.260971e+04 -4.014324e+03
56 -7.731240e+03 -1.260971e+04
57 -7.315045e+03 -7.731240e+03
58 -6.820506e+03 -7.315045e+03
59 -9.848030e+03 -6.820506e+03
60 -7.627914e+03 -9.848030e+03
61 -8.836213e+03 -7.627914e+03
62 -8.572741e+03 -8.836213e+03
63 -1.430316e+04 -8.572741e+03
64 -7.792093e+03 -1.430316e+04
65 -1.111790e+04 -7.792093e+03
66 -7.647476e+03 -1.111790e+04
67 -1.298757e+04 -7.647476e+03
68 -5.583541e+03 -1.298757e+04
69 -6.613939e+03 -5.583541e+03
70 -4.942520e+03 -6.613939e+03
71 -9.294798e+02 -4.942520e+03
72 5.409749e+03 -9.294798e+02
73 7.568573e+03 5.409749e+03
74 7.945646e+03 7.568573e+03
75 9.585543e+03 7.945646e+03
> 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/7t8ue1293031357.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/8t8ue1293031357.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/9miuz1293031357.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/10miuz1293031357.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/11pisn1293031357.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/12b19t1293031357.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/13pso11293031357.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/14stnp1293031357.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/15eblv1293031357.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/16hck11293031357.tab")
+ }
>
> try(system("convert tmp/1xyen1293031357.ps tmp/1xyen1293031357.png",intern=TRUE))
character(0)
> try(system("convert tmp/2q8e81293031357.ps tmp/2q8e81293031357.png",intern=TRUE))
character(0)
> try(system("convert tmp/3q8e81293031357.ps tmp/3q8e81293031357.png",intern=TRUE))
character(0)
> try(system("convert tmp/4q8e81293031357.ps tmp/4q8e81293031357.png",intern=TRUE))
character(0)
> try(system("convert tmp/5jzdb1293031357.ps tmp/5jzdb1293031357.png",intern=TRUE))
character(0)
> try(system("convert tmp/6jzdb1293031357.ps tmp/6jzdb1293031357.png",intern=TRUE))
character(0)
> try(system("convert tmp/7t8ue1293031357.ps tmp/7t8ue1293031357.png",intern=TRUE))
character(0)
> try(system("convert tmp/8t8ue1293031357.ps tmp/8t8ue1293031357.png",intern=TRUE))
character(0)
> try(system("convert tmp/9miuz1293031357.ps tmp/9miuz1293031357.png",intern=TRUE))
character(0)
> try(system("convert tmp/10miuz1293031357.ps tmp/10miuz1293031357.png",intern=TRUE))
character(0)
>
>
> proc.time()
user system elapsed
3.090 1.900 4.975