R version 2.15.2 (2012-10-26) -- "Trick or Treat"
Copyright (C) 2012 The R Foundation for Statistical Computing
ISBN 3-900051-07-0
Platform: i686-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(18.2
+ ,2687
+ ,1870
+ ,1890
+ ,145.7
+ ,352.2
+ ,143.8
+ ,13271
+ ,9115
+ ,8190
+ ,-279.0
+ ,83.0
+ ,23.4
+ ,13621
+ ,4848
+ ,4572
+ ,485.0
+ ,898.9
+ ,1.1
+ ,3614
+ ,367
+ ,90
+ ,14.1
+ ,24.6
+ ,49.5
+ ,6425
+ ,6131
+ ,2448
+ ,345.8
+ ,682.5
+ ,4.8
+ ,1022
+ ,1754
+ ,1370
+ ,72.0
+ ,119.5
+ ,20.8
+ ,1093
+ ,1679
+ ,1070
+ ,100.9
+ ,164.5
+ ,19.4
+ ,1529
+ ,1295
+ ,444
+ ,25.6
+ ,137.0
+ ,2.1
+ ,2788
+ ,271
+ ,304
+ ,23.5
+ ,28.9
+ ,79.4
+ ,19788
+ ,9084
+ ,10636
+ ,1092.9
+ ,2576.8
+ ,2.8
+ ,327
+ ,542
+ ,959
+ ,54.1
+ ,72.5
+ ,3.8
+ ,1117
+ ,1038
+ ,478
+ ,59.7
+ ,91.7
+ ,4.1
+ ,5401
+ ,550
+ ,376
+ ,25.6
+ ,37.5
+ ,13.2
+ ,1128
+ ,1516
+ ,430
+ ,-47.0
+ ,26.7
+ ,2.8
+ ,1633
+ ,701
+ ,679
+ ,74.3
+ ,135.9
+ ,48.5
+ ,44736
+ ,16197
+ ,4653
+ ,-732.5
+ ,-651.9
+ ,6.2
+ ,5651
+ ,1254
+ ,2002
+ ,310.7
+ ,407.9
+ ,10.8
+ ,5835
+ ,4053
+ ,1601
+ ,-93.8
+ ,173.8
+ ,3.8
+ ,278
+ ,205
+ ,853
+ ,44.8
+ ,50.5
+ ,21.9
+ ,5074
+ ,2557
+ ,1892
+ ,239.9
+ ,578.3
+ ,12.6
+ ,866
+ ,1487
+ ,944
+ ,71.7
+ ,115.4
+ ,128.0
+ ,4418
+ ,8793
+ ,4459
+ ,283.6
+ ,456.5
+ ,87.3
+ ,6914
+ ,7029
+ ,7957
+ ,400.6
+ ,754.7
+ ,16.0
+ ,862
+ ,1601
+ ,1093
+ ,66.9
+ ,106.8
+ ,0.7
+ ,401
+ ,176
+ ,1084
+ ,55.6
+ ,57.0
+ ,22.5
+ ,430
+ ,1155
+ ,1045
+ ,55.7
+ ,70.8
+ ,15.4
+ ,799
+ ,1140
+ ,683
+ ,57.6
+ ,89.2
+ ,3.0
+ ,4789
+ ,453
+ ,367
+ ,40.2
+ ,51.4
+ ,2.1
+ ,2548
+ ,264
+ ,181
+ ,22.2
+ ,26.2
+ ,4.1
+ ,5249
+ ,527
+ ,346
+ ,37.8
+ ,56.2
+ ,6.4
+ ,3494
+ ,1653
+ ,1442
+ ,160.9
+ ,320.3
+ ,26.6
+ ,1804
+ ,2564
+ ,483
+ ,70.5
+ ,164.9
+ ,304.0
+ ,26432
+ ,28285
+ ,33172
+ ,2336.0
+ ,3562.0
+ ,18.6
+ ,623
+ ,2247
+ ,797
+ ,57.0
+ ,93.8
+ ,65.0
+ ,1608
+ ,6615
+ ,829
+ ,56.1
+ ,134.0
+ ,66.2
+ ,4662
+ ,4781
+ ,2988
+ ,28.7
+ ,371.5
+ ,83.0
+ ,5769
+ ,6571
+ ,9462
+ ,482.0
+ ,792.0
+ ,62.0
+ ,6259
+ ,4152
+ ,3090
+ ,283.7
+ ,524.5
+ ,1.6
+ ,1654
+ ,451
+ ,779
+ ,84.8
+ ,130.4
+ ,400.2
+ ,52634
+ ,50056
+ ,95697
+ ,6555.0
+ ,9874.0
+ ,23.3
+ ,999
+ ,1878
+ ,393
+ ,-173.5
+ ,-108.1
+ ,4.6
+ ,1679
+ ,1354
+ ,687
+ ,93.8
+ ,154.6
+ ,164.6
+ ,4178
+ ,17124
+ ,2091
+ ,180.8
+ ,390.4
+ ,1.9
+ ,223
+ ,557
+ ,1040
+ ,60.6
+ ,63.7
+ ,57.5
+ ,6307
+ ,8199
+ ,598
+ ,-771.5
+ ,-524.3
+ ,2.4
+ ,3720
+ ,356
+ ,211
+ ,26.6
+ ,34.8
+ ,77.3
+ ,3442
+ ,5080
+ ,2673
+ ,235.4
+ ,361.5
+ ,15.8
+ ,33406
+ ,3222
+ ,1413
+ ,201.7
+ ,246.7
+ ,0.6
+ ,1257
+ ,355
+ ,181
+ ,167.5
+ ,304.0
+ ,3.5
+ ,1743
+ ,597
+ ,717
+ ,121.6
+ ,172.4
+ ,9.0
+ ,12505
+ ,1302
+ ,702
+ ,108.4
+ ,131.4
+ ,62.0
+ ,3940
+ ,4317
+ ,3940
+ ,315.2
+ ,566.3
+ ,7.4
+ ,8998
+ ,882
+ ,988
+ ,93.0
+ ,119.0
+ ,15.6
+ ,21419
+ ,2516
+ ,930
+ ,107.6
+ ,164.7
+ ,25.2
+ ,2366
+ ,3305
+ ,1117
+ ,131.2
+ ,256.5
+ ,25.4
+ ,2448
+ ,3484
+ ,1036
+ ,48.8
+ ,257.1
+ ,3.5
+ ,1440
+ ,1617
+ ,639
+ ,81.7
+ ,126.4
+ ,27.3
+ ,14045
+ ,15636
+ ,2754
+ ,418.0
+ ,1462.0
+ ,37.5
+ ,4084
+ ,4346
+ ,3023
+ ,302.7
+ ,521.7
+ ,3.4
+ ,3010
+ ,749
+ ,1120
+ ,146.3
+ ,209.2
+ ,14.3
+ ,1286
+ ,1734
+ ,361
+ ,69.2
+ ,145.7
+ ,6.1
+ ,707
+ ,706
+ ,275
+ ,61.4
+ ,77.8
+ ,4.9
+ ,3086
+ ,1739
+ ,1507
+ ,202.7
+ ,335.2
+ ,3.3
+ ,252
+ ,312
+ ,883
+ ,41.7
+ ,60.6
+ ,7.0
+ ,11052
+ ,1097
+ ,606
+ ,64.9
+ ,97.6
+ ,8.2
+ ,9672
+ ,1037
+ ,829
+ ,92.6
+ ,118.2
+ ,43.5
+ ,1112
+ ,3689
+ ,542
+ ,30.3
+ ,96.9
+ ,48.5
+ ,1104
+ ,5123
+ ,910
+ ,63.7
+ ,133.3
+ ,5.4
+ ,478
+ ,672
+ ,866
+ ,67.1
+ ,101.6
+ ,49.5
+ ,10348
+ ,5721
+ ,1915
+ ,223.6
+ ,322.5
+ ,29.1
+ ,2769
+ ,3725
+ ,663
+ ,-208.4
+ ,12.4
+ ,2.6
+ ,752
+ ,2149
+ ,101
+ ,11.1
+ ,15.2
+ ,0.8
+ ,4989
+ ,518
+ ,53
+ ,-3.1
+ ,-0.3
+ ,184.8
+ ,10528
+ ,14992
+ ,5377
+ ,312.7
+ ,710.7
+ ,2.3
+ ,1995
+ ,2662
+ ,341
+ ,34.7
+ ,100.7
+ ,8.0
+ ,2286
+ ,2235
+ ,2306
+ ,195.3
+ ,219.0
+ ,10.3
+ ,952
+ ,1307
+ ,309
+ ,35.4
+ ,92.8
+ ,50.0
+ ,2957
+ ,2806
+ ,457
+ ,40.6
+ ,93.5
+ ,118.1
+ ,2535
+ ,5958
+ ,1921
+ ,177.0
+ ,288.0)
+ ,dim=c(6
+ ,79)
+ ,dimnames=list(c('Aantal_werknemers'
+ ,'Activa'
+ ,'Omzet'
+ ,'Marktwaarde'
+ ,'Winst'
+ ,'Cashflow')
+ ,1:79))
> y <- array(NA,dim=c(6,79),dimnames=list(c('Aantal_werknemers','Activa','Omzet','Marktwaarde','Winst','Cashflow'),1:79))
> for (i in 1:dim(x)[1])
+ {
+ for (j in 1:dim(x)[2])
+ {
+ y[i,j] <- as.numeric(x[i,j])
+ }
+ }
> par20 = ''
> par19 = ''
> par18 = ''
> par17 = ''
> par16 = ''
> par15 = ''
> par14 = ''
> par13 = ''
> par12 = ''
> par11 = ''
> par10 = ''
> par9 = ''
> par8 = ''
> par7 = ''
> par6 = ''
> par5 = ''
> par4 = ''
> par3 = 'No Linear Trend'
> par2 = 'Do not include Seasonal Dummies'
> par1 = '1'
> library(lattice)
> library(lmtest)
Loading required package: zoo
Attaching package: 'zoo'
The following object(s) are masked from 'package:base':
as.Date, as.Date.numeric
> n25 <- 25 #minimum number of obs. for Goldfeld-Quandt test
> par1 <- as.numeric(par1)
> x <- t(y)
> k <- length(x[1,])
> n <- length(x[,1])
> x1 <- cbind(x[,par1], x[,1:k!=par1])
> mycolnames <- c(colnames(x)[par1], colnames(x)[1:k!=par1])
> colnames(x1) <- mycolnames #colnames(x)[par1]
> x <- x1
> if (par3 == 'First Differences'){
+ x2 <- array(0, dim=c(n-1,k), dimnames=list(1:(n-1), paste('(1-B)',colnames(x),sep='')))
+ for (i in 1:n-1) {
+ for (j in 1:k) {
+ x2[i,j] <- x[i+1,j] - x[i,j]
+ }
+ }
+ x <- x2
+ }
> if (par2 == 'Include Monthly Dummies'){
+ x2 <- array(0, dim=c(n,11), dimnames=list(1:n, paste('M', seq(1:11), sep ='')))
+ for (i in 1:11){
+ x2[seq(i,n,12),i] <- 1
+ }
+ x <- cbind(x, x2)
+ }
> if (par2 == 'Include Quarterly Dummies'){
+ x2 <- array(0, dim=c(n,3), dimnames=list(1:n, paste('Q', seq(1:3), sep ='')))
+ for (i in 1:3){
+ x2[seq(i,n,4),i] <- 1
+ }
+ x <- cbind(x, x2)
+ }
> k <- length(x[1,])
> if (par3 == 'Linear Trend'){
+ x <- cbind(x, c(1:n))
+ colnames(x)[k+1] <- 't'
+ }
> x
Aantal_werknemers Activa Omzet Marktwaarde Winst Cashflow
1 18.2 2687 1870 1890 145.7 352.2
2 143.8 13271 9115 8190 -279.0 83.0
3 23.4 13621 4848 4572 485.0 898.9
4 1.1 3614 367 90 14.1 24.6
5 49.5 6425 6131 2448 345.8 682.5
6 4.8 1022 1754 1370 72.0 119.5
7 20.8 1093 1679 1070 100.9 164.5
8 19.4 1529 1295 444 25.6 137.0
9 2.1 2788 271 304 23.5 28.9
10 79.4 19788 9084 10636 1092.9 2576.8
11 2.8 327 542 959 54.1 72.5
12 3.8 1117 1038 478 59.7 91.7
13 4.1 5401 550 376 25.6 37.5
14 13.2 1128 1516 430 -47.0 26.7
15 2.8 1633 701 679 74.3 135.9
16 48.5 44736 16197 4653 -732.5 -651.9
17 6.2 5651 1254 2002 310.7 407.9
18 10.8 5835 4053 1601 -93.8 173.8
19 3.8 278 205 853 44.8 50.5
20 21.9 5074 2557 1892 239.9 578.3
21 12.6 866 1487 944 71.7 115.4
22 128.0 4418 8793 4459 283.6 456.5
23 87.3 6914 7029 7957 400.6 754.7
24 16.0 862 1601 1093 66.9 106.8
25 0.7 401 176 1084 55.6 57.0
26 22.5 430 1155 1045 55.7 70.8
27 15.4 799 1140 683 57.6 89.2
28 3.0 4789 453 367 40.2 51.4
29 2.1 2548 264 181 22.2 26.2
30 4.1 5249 527 346 37.8 56.2
31 6.4 3494 1653 1442 160.9 320.3
32 26.6 1804 2564 483 70.5 164.9
33 304.0 26432 28285 33172 2336.0 3562.0
34 18.6 623 2247 797 57.0 93.8
35 65.0 1608 6615 829 56.1 134.0
36 66.2 4662 4781 2988 28.7 371.5
37 83.0 5769 6571 9462 482.0 792.0
38 62.0 6259 4152 3090 283.7 524.5
39 1.6 1654 451 779 84.8 130.4
40 400.2 52634 50056 95697 6555.0 9874.0
41 23.3 999 1878 393 -173.5 -108.1
42 4.6 1679 1354 687 93.8 154.6
43 164.6 4178 17124 2091 180.8 390.4
44 1.9 223 557 1040 60.6 63.7
45 57.5 6307 8199 598 -771.5 -524.3
46 2.4 3720 356 211 26.6 34.8
47 77.3 3442 5080 2673 235.4 361.5
48 15.8 33406 3222 1413 201.7 246.7
49 0.6 1257 355 181 167.5 304.0
50 3.5 1743 597 717 121.6 172.4
51 9.0 12505 1302 702 108.4 131.4
52 62.0 3940 4317 3940 315.2 566.3
53 7.4 8998 882 988 93.0 119.0
54 15.6 21419 2516 930 107.6 164.7
55 25.2 2366 3305 1117 131.2 256.5
56 25.4 2448 3484 1036 48.8 257.1
57 3.5 1440 1617 639 81.7 126.4
58 27.3 14045 15636 2754 418.0 1462.0
59 37.5 4084 4346 3023 302.7 521.7
60 3.4 3010 749 1120 146.3 209.2
61 14.3 1286 1734 361 69.2 145.7
62 6.1 707 706 275 61.4 77.8
63 4.9 3086 1739 1507 202.7 335.2
64 3.3 252 312 883 41.7 60.6
65 7.0 11052 1097 606 64.9 97.6
66 8.2 9672 1037 829 92.6 118.2
67 43.5 1112 3689 542 30.3 96.9
68 48.5 1104 5123 910 63.7 133.3
69 5.4 478 672 866 67.1 101.6
70 49.5 10348 5721 1915 223.6 322.5
71 29.1 2769 3725 663 -208.4 12.4
72 2.6 752 2149 101 11.1 15.2
73 0.8 4989 518 53 -3.1 -0.3
74 184.8 10528 14992 5377 312.7 710.7
75 2.3 1995 2662 341 34.7 100.7
76 8.0 2286 2235 2306 195.3 219.0
77 10.3 952 1307 309 35.4 92.8
78 50.0 2957 2806 457 40.6 93.5
79 118.1 2535 5958 1921 177.0 288.0
> k <- length(x[1,])
> df <- as.data.frame(x)
> (mylm <- lm(df))
Call:
lm(formula = df)
Coefficients:
(Intercept) Activa Omzet Marktwaarde Winst Cashflow
6.4546652 -0.0014858 0.0104250 0.0006733 0.0445309 -0.0377493
> (mysum <- summary(mylm))
Call:
lm(formula = df)
Residuals:
Min 1Q Median 3Q Max
-86.570 -8.093 -3.793 4.946 72.083
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 6.4546652 3.4430353 1.875 0.06483 .
Activa -0.0014858 0.0004391 -3.384 0.00115 **
Omzet 0.0104250 0.0009730 10.715 < 2e-16 ***
Marktwaarde 0.0006733 0.0012076 0.558 0.57886
Winst 0.0445309 0.0276867 1.608 0.11207
Cashflow -0.0377493 0.0177168 -2.131 0.03648 *
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 22.84 on 73 degrees of freedom
Multiple R-squared: 0.8827, Adjusted R-squared: 0.8746
F-statistic: 109.8 on 5 and 73 DF, p-value: < 2.2e-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.0534622903 0.1069245807 0.94653771
[2,] 0.0222399831 0.0444799662 0.97776002
[3,] 0.0064415271 0.0128830542 0.99355847
[4,] 0.0019887262 0.0039774523 0.99801127
[5,] 0.0005012808 0.0010025616 0.99949872
[6,] 0.0008407663 0.0016815326 0.99915923
[7,] 0.0002378133 0.0004756266 0.99976219
[8,] 0.0426728168 0.0853456336 0.95732718
[9,] 0.0252963918 0.0505927836 0.97470361
[10,] 0.0434448634 0.0868897268 0.95655514
[11,] 0.0283401265 0.0566802529 0.97165987
[12,] 0.0260229632 0.0520459263 0.97397704
[13,] 0.0156791767 0.0313583534 0.98432082
[14,] 0.0232579253 0.0465158505 0.97674207
[15,] 0.0949308454 0.1898616908 0.90506915
[16,] 0.0668486593 0.1336973186 0.93315134
[17,] 0.0485656269 0.0971312539 0.95143437
[18,] 0.0325427252 0.0650854505 0.96745727
[19,] 0.0206330823 0.0412661646 0.97936692
[20,] 0.0190063286 0.0380126571 0.98099367
[21,] 0.0130012638 0.0260025276 0.98699874
[22,] 0.0119134773 0.0238269546 0.98808652
[23,] 0.0076683011 0.0153366021 0.99233170
[24,] 0.0047351336 0.0094702672 0.99526487
[25,] 0.0283651044 0.0567302088 0.97163490
[26,] 0.0215012577 0.0430025154 0.97849874
[27,] 0.0147414686 0.0294829372 0.98525853
[28,] 0.0254504843 0.0509009686 0.97454952
[29,] 0.0212637351 0.0425274702 0.97873626
[30,] 0.0679045228 0.1358090456 0.93209548
[31,] 0.0482335857 0.0964671713 0.95176641
[32,] 0.7421006560 0.5157986880 0.25789934
[33,] 0.7052505623 0.5894988754 0.29474944
[34,] 0.6572903591 0.6854192818 0.34270964
[35,] 0.8455016317 0.3089967365 0.15449837
[36,] 0.8217612201 0.3564775597 0.17823878
[37,] 0.9061766641 0.1876466718 0.09382334
[38,] 0.8735460391 0.2529079219 0.12645396
[39,] 0.8652314540 0.2695370920 0.13476855
[40,] 0.9045201027 0.1909597946 0.09547990
[41,] 0.9346422391 0.1307155217 0.06535776
[42,] 0.9111179671 0.1777640658 0.08888203
[43,] 0.8815050528 0.2369898943 0.11849495
[44,] 0.8754049473 0.2491901053 0.12459505
[45,] 0.8334233589 0.3331532821 0.16657664
[46,] 0.7916476084 0.4167047832 0.20835239
[47,] 0.7364783956 0.5270432089 0.26352160
[48,] 0.6874394305 0.6251211391 0.31256057
[49,] 0.6379348362 0.7241303277 0.36206516
[50,] 0.9090021708 0.1819956584 0.09099783
[51,] 0.8727660877 0.2544678245 0.12723391
[52,] 0.8195601676 0.3608796647 0.18043983
[53,] 0.7515068488 0.4969863024 0.24849315
[54,] 0.6730167383 0.6539665235 0.32698326
[55,] 0.6299782095 0.7400435809 0.37002179
[56,] 0.5398092355 0.9203815290 0.46019076
[57,] 0.4450492572 0.8900985144 0.55495074
[58,] 0.3755796919 0.7511593838 0.62442031
[59,] 0.2700496622 0.5400993245 0.72995034
[60,] 0.2087258819 0.4174517638 0.79127412
[61,] 0.1275977336 0.2551954672 0.87240227
[62,] 0.0768315760 0.1536631520 0.92316842
> postscript(file="/var/wessaorg/rcomp/tmp/1ahy01355676657.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/2q7ud1355676657.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/34n6y1355676657.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/485up1355676657.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/5qare1355676657.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 = 79
Frequency = 1
1 2 3 4 5 6
1.77768347 72.08314968 -4.09932016 -3.57065506 -2.60681418 -18.03916470
7 8 9 10 11 12
-0.53802642 5.44955644 -3.19751857 49.08996515 -9.13713492 -11.33483298
13 14 15 16 17 18
-0.04093545 -4.57156689 -7.17188781 -55.46029176 -4.71687426 -19.57735730
19 20 21 22 23 24
-5.04170043 6.20139170 -7.54206894 38.04425431 23.13406565 -5.54767916
25 26 27 28 29 30
-8.04771594 4.13206691 -1.40956530 -1.15842273 -3.44236112 0.15582611
31 32 33 34 35 36
-8.14045365 -1.14362566 50.05239401 -9.88791594 -6.02461187 27.56449169
37 38 39 40 41 42
18.67731749 26.64627585 -6.47696793 -33.47873408 2.13238060 -12.27887664
43 44 45 46 47 48
-8.88589037 -11.02422050 -10.89687730 -2.25154518 24.36475513 24.77152684
49 50 51 52 53 54
-3.79282959 -5.97829520 7.21289461 21.08337409 4.80566031 15.54091562
55 56 57 58 59 60
-9.10555888 -6.90325485 -16.96915011 -86.56967844 -4.01454671 -5.76241462
61 62 63 64 65 66
-6.14531239 -6.64667104 -12.48591740 -6.19669003 5.91689433 5.08591459
67 68 69 70 71 72
2.18353480 -8.13880596 -7.08579828 -0.29283120 -2.77148507 -25.12909880
73 74 75 76 77 78
-3.55091041 46.98023367 -26.91516854 -20.34029993 -6.64690030 20.10038120
79
54.99632667
> postscript(file="/var/wessaorg/rcomp/tmp/6yrmn1355676657.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 = 79
Frequency = 1
lag(myerror, k = 1) myerror
0 1.77768347 NA
1 72.08314968 1.77768347
2 -4.09932016 72.08314968
3 -3.57065506 -4.09932016
4 -2.60681418 -3.57065506
5 -18.03916470 -2.60681418
6 -0.53802642 -18.03916470
7 5.44955644 -0.53802642
8 -3.19751857 5.44955644
9 49.08996515 -3.19751857
10 -9.13713492 49.08996515
11 -11.33483298 -9.13713492
12 -0.04093545 -11.33483298
13 -4.57156689 -0.04093545
14 -7.17188781 -4.57156689
15 -55.46029176 -7.17188781
16 -4.71687426 -55.46029176
17 -19.57735730 -4.71687426
18 -5.04170043 -19.57735730
19 6.20139170 -5.04170043
20 -7.54206894 6.20139170
21 38.04425431 -7.54206894
22 23.13406565 38.04425431
23 -5.54767916 23.13406565
24 -8.04771594 -5.54767916
25 4.13206691 -8.04771594
26 -1.40956530 4.13206691
27 -1.15842273 -1.40956530
28 -3.44236112 -1.15842273
29 0.15582611 -3.44236112
30 -8.14045365 0.15582611
31 -1.14362566 -8.14045365
32 50.05239401 -1.14362566
33 -9.88791594 50.05239401
34 -6.02461187 -9.88791594
35 27.56449169 -6.02461187
36 18.67731749 27.56449169
37 26.64627585 18.67731749
38 -6.47696793 26.64627585
39 -33.47873408 -6.47696793
40 2.13238060 -33.47873408
41 -12.27887664 2.13238060
42 -8.88589037 -12.27887664
43 -11.02422050 -8.88589037
44 -10.89687730 -11.02422050
45 -2.25154518 -10.89687730
46 24.36475513 -2.25154518
47 24.77152684 24.36475513
48 -3.79282959 24.77152684
49 -5.97829520 -3.79282959
50 7.21289461 -5.97829520
51 21.08337409 7.21289461
52 4.80566031 21.08337409
53 15.54091562 4.80566031
54 -9.10555888 15.54091562
55 -6.90325485 -9.10555888
56 -16.96915011 -6.90325485
57 -86.56967844 -16.96915011
58 -4.01454671 -86.56967844
59 -5.76241462 -4.01454671
60 -6.14531239 -5.76241462
61 -6.64667104 -6.14531239
62 -12.48591740 -6.64667104
63 -6.19669003 -12.48591740
64 5.91689433 -6.19669003
65 5.08591459 5.91689433
66 2.18353480 5.08591459
67 -8.13880596 2.18353480
68 -7.08579828 -8.13880596
69 -0.29283120 -7.08579828
70 -2.77148507 -0.29283120
71 -25.12909880 -2.77148507
72 -3.55091041 -25.12909880
73 46.98023367 -3.55091041
74 -26.91516854 46.98023367
75 -20.34029993 -26.91516854
76 -6.64690030 -20.34029993
77 20.10038120 -6.64690030
78 54.99632667 20.10038120
79 NA 54.99632667
> dum1 <- dum[2:length(myerror),]
> dum1
lag(myerror, k = 1) myerror
[1,] 72.08314968 1.77768347
[2,] -4.09932016 72.08314968
[3,] -3.57065506 -4.09932016
[4,] -2.60681418 -3.57065506
[5,] -18.03916470 -2.60681418
[6,] -0.53802642 -18.03916470
[7,] 5.44955644 -0.53802642
[8,] -3.19751857 5.44955644
[9,] 49.08996515 -3.19751857
[10,] -9.13713492 49.08996515
[11,] -11.33483298 -9.13713492
[12,] -0.04093545 -11.33483298
[13,] -4.57156689 -0.04093545
[14,] -7.17188781 -4.57156689
[15,] -55.46029176 -7.17188781
[16,] -4.71687426 -55.46029176
[17,] -19.57735730 -4.71687426
[18,] -5.04170043 -19.57735730
[19,] 6.20139170 -5.04170043
[20,] -7.54206894 6.20139170
[21,] 38.04425431 -7.54206894
[22,] 23.13406565 38.04425431
[23,] -5.54767916 23.13406565
[24,] -8.04771594 -5.54767916
[25,] 4.13206691 -8.04771594
[26,] -1.40956530 4.13206691
[27,] -1.15842273 -1.40956530
[28,] -3.44236112 -1.15842273
[29,] 0.15582611 -3.44236112
[30,] -8.14045365 0.15582611
[31,] -1.14362566 -8.14045365
[32,] 50.05239401 -1.14362566
[33,] -9.88791594 50.05239401
[34,] -6.02461187 -9.88791594
[35,] 27.56449169 -6.02461187
[36,] 18.67731749 27.56449169
[37,] 26.64627585 18.67731749
[38,] -6.47696793 26.64627585
[39,] -33.47873408 -6.47696793
[40,] 2.13238060 -33.47873408
[41,] -12.27887664 2.13238060
[42,] -8.88589037 -12.27887664
[43,] -11.02422050 -8.88589037
[44,] -10.89687730 -11.02422050
[45,] -2.25154518 -10.89687730
[46,] 24.36475513 -2.25154518
[47,] 24.77152684 24.36475513
[48,] -3.79282959 24.77152684
[49,] -5.97829520 -3.79282959
[50,] 7.21289461 -5.97829520
[51,] 21.08337409 7.21289461
[52,] 4.80566031 21.08337409
[53,] 15.54091562 4.80566031
[54,] -9.10555888 15.54091562
[55,] -6.90325485 -9.10555888
[56,] -16.96915011 -6.90325485
[57,] -86.56967844 -16.96915011
[58,] -4.01454671 -86.56967844
[59,] -5.76241462 -4.01454671
[60,] -6.14531239 -5.76241462
[61,] -6.64667104 -6.14531239
[62,] -12.48591740 -6.64667104
[63,] -6.19669003 -12.48591740
[64,] 5.91689433 -6.19669003
[65,] 5.08591459 5.91689433
[66,] 2.18353480 5.08591459
[67,] -8.13880596 2.18353480
[68,] -7.08579828 -8.13880596
[69,] -0.29283120 -7.08579828
[70,] -2.77148507 -0.29283120
[71,] -25.12909880 -2.77148507
[72,] -3.55091041 -25.12909880
[73,] 46.98023367 -3.55091041
[74,] -26.91516854 46.98023367
[75,] -20.34029993 -26.91516854
[76,] -6.64690030 -20.34029993
[77,] 20.10038120 -6.64690030
[78,] 54.99632667 20.10038120
> z <- as.data.frame(dum1)
> z
lag(myerror, k = 1) myerror
1 72.08314968 1.77768347
2 -4.09932016 72.08314968
3 -3.57065506 -4.09932016
4 -2.60681418 -3.57065506
5 -18.03916470 -2.60681418
6 -0.53802642 -18.03916470
7 5.44955644 -0.53802642
8 -3.19751857 5.44955644
9 49.08996515 -3.19751857
10 -9.13713492 49.08996515
11 -11.33483298 -9.13713492
12 -0.04093545 -11.33483298
13 -4.57156689 -0.04093545
14 -7.17188781 -4.57156689
15 -55.46029176 -7.17188781
16 -4.71687426 -55.46029176
17 -19.57735730 -4.71687426
18 -5.04170043 -19.57735730
19 6.20139170 -5.04170043
20 -7.54206894 6.20139170
21 38.04425431 -7.54206894
22 23.13406565 38.04425431
23 -5.54767916 23.13406565
24 -8.04771594 -5.54767916
25 4.13206691 -8.04771594
26 -1.40956530 4.13206691
27 -1.15842273 -1.40956530
28 -3.44236112 -1.15842273
29 0.15582611 -3.44236112
30 -8.14045365 0.15582611
31 -1.14362566 -8.14045365
32 50.05239401 -1.14362566
33 -9.88791594 50.05239401
34 -6.02461187 -9.88791594
35 27.56449169 -6.02461187
36 18.67731749 27.56449169
37 26.64627585 18.67731749
38 -6.47696793 26.64627585
39 -33.47873408 -6.47696793
40 2.13238060 -33.47873408
41 -12.27887664 2.13238060
42 -8.88589037 -12.27887664
43 -11.02422050 -8.88589037
44 -10.89687730 -11.02422050
45 -2.25154518 -10.89687730
46 24.36475513 -2.25154518
47 24.77152684 24.36475513
48 -3.79282959 24.77152684
49 -5.97829520 -3.79282959
50 7.21289461 -5.97829520
51 21.08337409 7.21289461
52 4.80566031 21.08337409
53 15.54091562 4.80566031
54 -9.10555888 15.54091562
55 -6.90325485 -9.10555888
56 -16.96915011 -6.90325485
57 -86.56967844 -16.96915011
58 -4.01454671 -86.56967844
59 -5.76241462 -4.01454671
60 -6.14531239 -5.76241462
61 -6.64667104 -6.14531239
62 -12.48591740 -6.64667104
63 -6.19669003 -12.48591740
64 5.91689433 -6.19669003
65 5.08591459 5.91689433
66 2.18353480 5.08591459
67 -8.13880596 2.18353480
68 -7.08579828 -8.13880596
69 -0.29283120 -7.08579828
70 -2.77148507 -0.29283120
71 -25.12909880 -2.77148507
72 -3.55091041 -25.12909880
73 46.98023367 -3.55091041
74 -26.91516854 46.98023367
75 -20.34029993 -26.91516854
76 -6.64690030 -20.34029993
77 20.10038120 -6.64690030
78 54.99632667 20.10038120
> 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/73x931355676657.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/8tuu51355676657.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/966rt1355676657.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/10tksp1355676657.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/11i9eu1355676657.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/12c2dy1355676657.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/13cndo1355676657.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/14wkma1355676657.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/15qe731355676657.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/16w4v61355676657.tab")
+ }
>
> try(system("convert tmp/1ahy01355676657.ps tmp/1ahy01355676657.png",intern=TRUE))
character(0)
> try(system("convert tmp/2q7ud1355676657.ps tmp/2q7ud1355676657.png",intern=TRUE))
character(0)
> try(system("convert tmp/34n6y1355676657.ps tmp/34n6y1355676657.png",intern=TRUE))
character(0)
> try(system("convert tmp/485up1355676657.ps tmp/485up1355676657.png",intern=TRUE))
character(0)
> try(system("convert tmp/5qare1355676657.ps tmp/5qare1355676657.png",intern=TRUE))
character(0)
> try(system("convert tmp/6yrmn1355676657.ps tmp/6yrmn1355676657.png",intern=TRUE))
character(0)
> try(system("convert tmp/73x931355676657.ps tmp/73x931355676657.png",intern=TRUE))
character(0)
> try(system("convert tmp/8tuu51355676657.ps tmp/8tuu51355676657.png",intern=TRUE))
character(0)
> try(system("convert tmp/966rt1355676657.ps tmp/966rt1355676657.png",intern=TRUE))
character(0)
> try(system("convert tmp/10tksp1355676657.ps tmp/10tksp1355676657.png",intern=TRUE))
character(0)
>
>
> proc.time()
user system elapsed
6.143 1.192 7.368