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
+ ,143.8
+ ,13271
+ ,9115
+ ,8190
+ ,23.4
+ ,13621
+ ,4848
+ ,4572
+ ,1.1
+ ,3614
+ ,367
+ ,90
+ ,49.5
+ ,6425
+ ,6131
+ ,2448
+ ,4.8
+ ,1022
+ ,1754
+ ,1370
+ ,20.8
+ ,1093
+ ,1679
+ ,1070
+ ,19.4
+ ,1529
+ ,1295
+ ,444
+ ,2.1
+ ,2788
+ ,271
+ ,304
+ ,79.4
+ ,19788
+ ,9084
+ ,10636
+ ,2.8
+ ,327
+ ,542
+ ,959
+ ,3.8
+ ,1117
+ ,1038
+ ,478
+ ,4.1
+ ,5401
+ ,550
+ ,376
+ ,13.2
+ ,1128
+ ,1516
+ ,430
+ ,2.8
+ ,1633
+ ,701
+ ,679
+ ,48.5
+ ,44736
+ ,16197
+ ,4653
+ ,6.2
+ ,5651
+ ,1254
+ ,2002
+ ,10.8
+ ,5835
+ ,4053
+ ,1601
+ ,3.8
+ ,278
+ ,205
+ ,853
+ ,21.9
+ ,5074
+ ,2557
+ ,1892
+ ,12.6
+ ,866
+ ,1487
+ ,944
+ ,128.0
+ ,4418
+ ,8793
+ ,4459
+ ,87.3
+ ,6914
+ ,7029
+ ,7957
+ ,16.0
+ ,862
+ ,1601
+ ,1093
+ ,0.7
+ ,401
+ ,176
+ ,1084
+ ,22.5
+ ,430
+ ,1155
+ ,1045
+ ,15.4
+ ,799
+ ,1140
+ ,683
+ ,3.0
+ ,4789
+ ,453
+ ,367
+ ,2.1
+ ,2548
+ ,264
+ ,181
+ ,4.1
+ ,5249
+ ,527
+ ,346
+ ,6.4
+ ,3494
+ ,1653
+ ,1442
+ ,26.6
+ ,1804
+ ,2564
+ ,483
+ ,304.0
+ ,26432
+ ,28285
+ ,33172
+ ,18.6
+ ,623
+ ,2247
+ ,797
+ ,65.0
+ ,1608
+ ,6615
+ ,829
+ ,66.2
+ ,4662
+ ,4781
+ ,2988
+ ,83.0
+ ,5769
+ ,6571
+ ,9462
+ ,62.0
+ ,6259
+ ,4152
+ ,3090
+ ,1.6
+ ,1654
+ ,451
+ ,779
+ ,400.2
+ ,52634
+ ,50056
+ ,95697
+ ,23.3
+ ,999
+ ,1878
+ ,393
+ ,4.6
+ ,1679
+ ,1354
+ ,687
+ ,164.6
+ ,4178
+ ,17124
+ ,2091
+ ,1.9
+ ,223
+ ,557
+ ,1040
+ ,57.5
+ ,6307
+ ,8199
+ ,598
+ ,2.4
+ ,3720
+ ,356
+ ,211
+ ,77.3
+ ,3442
+ ,5080
+ ,2673
+ ,15.8
+ ,33406
+ ,3222
+ ,1413
+ ,0.6
+ ,1257
+ ,355
+ ,181
+ ,3.5
+ ,1743
+ ,597
+ ,717
+ ,9.0
+ ,12505
+ ,1302
+ ,702
+ ,62.0
+ ,3940
+ ,4317
+ ,3940
+ ,7.4
+ ,8998
+ ,882
+ ,988
+ ,15.6
+ ,21419
+ ,2516
+ ,930
+ ,25.2
+ ,2366
+ ,3305
+ ,1117
+ ,25.4
+ ,2448
+ ,3484
+ ,1036
+ ,3.5
+ ,1440
+ ,1617
+ ,639
+ ,27.3
+ ,14045
+ ,15636
+ ,2754
+ ,37.5
+ ,4084
+ ,4346
+ ,3023
+ ,3.4
+ ,3010
+ ,749
+ ,1120
+ ,14.3
+ ,1286
+ ,1734
+ ,361
+ ,6.1
+ ,707
+ ,706
+ ,275
+ ,4.9
+ ,3086
+ ,1739
+ ,1507
+ ,3.3
+ ,252
+ ,312
+ ,883
+ ,7.0
+ ,11052
+ ,1097
+ ,606
+ ,8.2
+ ,9672
+ ,1037
+ ,829
+ ,43.5
+ ,1112
+ ,3689
+ ,542
+ ,48.5
+ ,1104
+ ,5123
+ ,910
+ ,5.4
+ ,478
+ ,672
+ ,866
+ ,49.5
+ ,10348
+ ,5721
+ ,1915
+ ,29.1
+ ,2769
+ ,3725
+ ,663
+ ,2.6
+ ,752
+ ,2149
+ ,101
+ ,0.8
+ ,4989
+ ,518
+ ,53
+ ,184.8
+ ,10528
+ ,14992
+ ,5377
+ ,2.3
+ ,1995
+ ,2662
+ ,341
+ ,8.0
+ ,2286
+ ,2235
+ ,2306
+ ,10.3
+ ,952
+ ,1307
+ ,309
+ ,50.0
+ ,2957
+ ,2806
+ ,457
+ ,118.1
+ ,2535
+ ,5958
+ ,1921)
+ ,dim=c(4
+ ,79)
+ ,dimnames=list(c('Aantal_werknemers'
+ ,'Activa'
+ ,'Omzet'
+ ,'Marktwaarde')
+ ,1:79))
> y <- array(NA,dim=c(4,79),dimnames=list(c('Aantal_werknemers','Activa','Omzet','Marktwaarde'),1:79))
> 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'
> 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
1 18.2 2687 1870 1890
2 143.8 13271 9115 8190
3 23.4 13621 4848 4572
4 1.1 3614 367 90
5 49.5 6425 6131 2448
6 4.8 1022 1754 1370
7 20.8 1093 1679 1070
8 19.4 1529 1295 444
9 2.1 2788 271 304
10 79.4 19788 9084 10636
11 2.8 327 542 959
12 3.8 1117 1038 478
13 4.1 5401 550 376
14 13.2 1128 1516 430
15 2.8 1633 701 679
16 48.5 44736 16197 4653
17 6.2 5651 1254 2002
18 10.8 5835 4053 1601
19 3.8 278 205 853
20 21.9 5074 2557 1892
21 12.6 866 1487 944
22 128.0 4418 8793 4459
23 87.3 6914 7029 7957
24 16.0 862 1601 1093
25 0.7 401 176 1084
26 22.5 430 1155 1045
27 15.4 799 1140 683
28 3.0 4789 453 367
29 2.1 2548 264 181
30 4.1 5249 527 346
31 6.4 3494 1653 1442
32 26.6 1804 2564 483
33 304.0 26432 28285 33172
34 18.6 623 2247 797
35 65.0 1608 6615 829
36 66.2 4662 4781 2988
37 83.0 5769 6571 9462
38 62.0 6259 4152 3090
39 1.6 1654 451 779
40 400.2 52634 50056 95697
41 23.3 999 1878 393
42 4.6 1679 1354 687
43 164.6 4178 17124 2091
44 1.9 223 557 1040
45 57.5 6307 8199 598
46 2.4 3720 356 211
47 77.3 3442 5080 2673
48 15.8 33406 3222 1413
49 0.6 1257 355 181
50 3.5 1743 597 717
51 9.0 12505 1302 702
52 62.0 3940 4317 3940
53 7.4 8998 882 988
54 15.6 21419 2516 930
55 25.2 2366 3305 1117
56 25.4 2448 3484 1036
57 3.5 1440 1617 639
58 27.3 14045 15636 2754
59 37.5 4084 4346 3023
60 3.4 3010 749 1120
61 14.3 1286 1734 361
62 6.1 707 706 275
63 4.9 3086 1739 1507
64 3.3 252 312 883
65 7.0 11052 1097 606
66 8.2 9672 1037 829
67 43.5 1112 3689 542
68 48.5 1104 5123 910
69 5.4 478 672 866
70 49.5 10348 5721 1915
71 29.1 2769 3725 663
72 2.6 752 2149 101
73 0.8 4989 518 53
74 184.8 10528 14992 5377
75 2.3 1995 2662 341
76 8.0 2286 2235 2306
77 10.3 952 1307 309
78 50.0 2957 2806 457
79 118.1 2535 5958 1921
> k <- length(x[1,])
> df <- as.data.frame(x)
> (mylm <- lm(df))
Call:
lm(formula = df)
Coefficients:
(Intercept) Activa Omzet Marktwaarde
5.7373179 -0.0015394 0.0095823 0.0002958
> (mysum <- summary(mylm))
Call:
lm(formula = df)
Residuals:
Min 1Q Median 3Q Max
-107.460 -7.928 -3.447 5.177 68.726
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 5.7373179 3.4385321 1.669 0.099379 .
Activa -0.0015394 0.0004336 -3.550 0.000669 ***
Omzet 0.0095823 0.0008680 11.040 < 2e-16 ***
Marktwaarde 0.0002958 0.0004900 0.604 0.547829
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 23.26 on 75 degrees of freedom
Multiple R-squared: 0.875, Adjusted R-squared: 0.87
F-statistic: 175 on 3 and 75 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.4919609576 0.9839219152 0.50803904
[2,] 0.4014270606 0.8028541213 0.59857294
[3,] 0.3015164157 0.6030328314 0.69848358
[4,] 0.4073884202 0.8147768404 0.59261158
[5,] 0.3022473951 0.6044947902 0.69775260
[6,] 0.2119994291 0.4239988581 0.78800057
[7,] 0.1877509598 0.3755019197 0.81224904
[8,] 0.1239170623 0.2478341245 0.87608294
[9,] 0.0783958108 0.1567916216 0.92160419
[10,] 0.0662969830 0.1325939661 0.93370302
[11,] 0.0401576122 0.0803152244 0.95984239
[12,] 0.0497126500 0.0994252999 0.95028735
[13,] 0.0301284494 0.0602568988 0.96987155
[14,] 0.0176735959 0.0353471918 0.98232640
[15,] 0.0102911352 0.0205822704 0.98970886
[16,] 0.0116323492 0.0232646984 0.98836765
[17,] 0.0115956714 0.0231913428 0.98840433
[18,] 0.0067648557 0.0135297113 0.99323514
[19,] 0.0038118588 0.0076237177 0.99618814
[20,] 0.0022561419 0.0045122838 0.99774386
[21,] 0.0012042089 0.0024084177 0.99879579
[22,] 0.0009763644 0.0019527288 0.99902364
[23,] 0.0005725135 0.0011450271 0.99942749
[24,] 0.0004549451 0.0009098901 0.99954505
[25,] 0.0002632862 0.0005265725 0.99973671
[26,] 0.0001285361 0.0002570723 0.99987146
[27,] 0.0046598423 0.0093196847 0.99534016
[28,] 0.0031159951 0.0062319902 0.99688400
[29,] 0.0019096305 0.0038192611 0.99809037
[30,] 0.0018657073 0.0037314146 0.99813429
[31,] 0.0013241532 0.0026483063 0.99867585
[32,] 0.0019031428 0.0038062857 0.99809686
[33,] 0.0011192034 0.0022384068 0.99888080
[34,] 0.0627599374 0.1255198747 0.93724006
[35,] 0.0462379373 0.0924758747 0.95376206
[36,] 0.0357880151 0.0715760303 0.96421198
[37,] 0.1325212629 0.2650425258 0.86747874
[38,] 0.1084421602 0.2168843204 0.89155784
[39,] 0.1133340163 0.2266680327 0.88666598
[40,] 0.0841359423 0.1682718847 0.91586406
[41,] 0.0812050280 0.1624100560 0.91879497
[42,] 0.1223059453 0.2446118907 0.87769405
[43,] 0.0918490505 0.1836981010 0.90815095
[44,] 0.0680672271 0.1361344542 0.93193277
[45,] 0.0512747606 0.1025495212 0.94872524
[46,] 0.0419865924 0.0839731847 0.95801341
[47,] 0.0290359673 0.0580719345 0.97096403
[48,] 0.0288908971 0.0577817943 0.97110910
[49,] 0.0201061581 0.0402123163 0.97989384
[50,] 0.0137601064 0.0275202127 0.98623989
[51,] 0.0102969884 0.0205939767 0.98970301
[52,] 0.9535937511 0.0928124977 0.04640625
[53,] 0.9361058105 0.1277883790 0.06389419
[54,] 0.9043183321 0.1913633359 0.09568167
[55,] 0.8608157556 0.2783684889 0.13918424
[56,] 0.8073197361 0.3853605278 0.19268026
[57,] 0.7571253651 0.4857492698 0.24287463
[58,] 0.6807618231 0.6384763537 0.31923818
[59,] 0.5962651994 0.8074696013 0.40373480
[60,] 0.5274728957 0.9450542087 0.47252710
[61,] 0.4226351415 0.8452702830 0.57736486
[62,] 0.3595559415 0.7191118830 0.64044406
[63,] 0.2604429143 0.5208858287 0.73955709
[64,] 0.1736444384 0.3472888768 0.82635556
[65,] 0.1183849303 0.2367698606 0.88161507
[66,] 0.1187466381 0.2374932761 0.88125336
> postscript(file="/var/wessaorg/rcomp/tmp/1f0dj1355997553.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/229d61355997553.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/3ollp1355997553.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/405l01355997553.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/5z6sw1355997553.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.8790100 68.7263705 -9.1766451 -2.6171986 -5.8200672 -16.5767531
7 8 9 10 11 12
0.3399743 3.4759765 -2.0321659 13.9321823 -7.9112684 -10.3056632
13 14 15 16 17 18
1.2955652 -5.4548847 -7.3415410 -44.9514814 -3.4465874 -25.2656536
19 20 21 22 23 24
-3.7260937 -1.0880709 -6.3323893 43.4871998 22.4979617 -4.0750140
25 26 27 28 29 30
-6.4271977 6.0478784 -0.2332449 0.1855895 -2.2981612 1.2908424
31 32 33 34 35 36
-10.2247933 -1.0722031 58.1025208 -7.9455529 -1.8943273 20.9423309
37 38 39 40 41 42
20.3787885 25.1978958 -6.1432142 -32.4763259 0.9886715 -11.7303586
43 44 45 46 47 48
0.5878749 -9.1390664 -17.2706721 -1.0844117 27.3923124 30.1962195
49 50 51 52 53 54
-6.6575442 -5.4868842 9.8292404 19.7954284 6.7704649 18.4512193
55 56 57 58 59 60
-8.8951248 -10.2601667 -15.7042332 -107.4603053 -4.4894928 -5.2121801
61 62 63 64 65 66
-6.1801919 -5.3954336 -13.1961871 -5.3003037 7.5852435 7.1698112
67 68 69 70 71 72
3.9649379 -4.8973156 -6.2970058 4.3055179 -8.2650053 -22.6019908
73 74 75 76 77 78
-2.2364830 50.0205777 -25.9752327 -16.3169407 -6.5873174 21.7915147
79
58.6052452
> postscript(file="/var/wessaorg/rcomp/tmp/65zc81355997553.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.8790100 NA
1 68.7263705 -1.8790100
2 -9.1766451 68.7263705
3 -2.6171986 -9.1766451
4 -5.8200672 -2.6171986
5 -16.5767531 -5.8200672
6 0.3399743 -16.5767531
7 3.4759765 0.3399743
8 -2.0321659 3.4759765
9 13.9321823 -2.0321659
10 -7.9112684 13.9321823
11 -10.3056632 -7.9112684
12 1.2955652 -10.3056632
13 -5.4548847 1.2955652
14 -7.3415410 -5.4548847
15 -44.9514814 -7.3415410
16 -3.4465874 -44.9514814
17 -25.2656536 -3.4465874
18 -3.7260937 -25.2656536
19 -1.0880709 -3.7260937
20 -6.3323893 -1.0880709
21 43.4871998 -6.3323893
22 22.4979617 43.4871998
23 -4.0750140 22.4979617
24 -6.4271977 -4.0750140
25 6.0478784 -6.4271977
26 -0.2332449 6.0478784
27 0.1855895 -0.2332449
28 -2.2981612 0.1855895
29 1.2908424 -2.2981612
30 -10.2247933 1.2908424
31 -1.0722031 -10.2247933
32 58.1025208 -1.0722031
33 -7.9455529 58.1025208
34 -1.8943273 -7.9455529
35 20.9423309 -1.8943273
36 20.3787885 20.9423309
37 25.1978958 20.3787885
38 -6.1432142 25.1978958
39 -32.4763259 -6.1432142
40 0.9886715 -32.4763259
41 -11.7303586 0.9886715
42 0.5878749 -11.7303586
43 -9.1390664 0.5878749
44 -17.2706721 -9.1390664
45 -1.0844117 -17.2706721
46 27.3923124 -1.0844117
47 30.1962195 27.3923124
48 -6.6575442 30.1962195
49 -5.4868842 -6.6575442
50 9.8292404 -5.4868842
51 19.7954284 9.8292404
52 6.7704649 19.7954284
53 18.4512193 6.7704649
54 -8.8951248 18.4512193
55 -10.2601667 -8.8951248
56 -15.7042332 -10.2601667
57 -107.4603053 -15.7042332
58 -4.4894928 -107.4603053
59 -5.2121801 -4.4894928
60 -6.1801919 -5.2121801
61 -5.3954336 -6.1801919
62 -13.1961871 -5.3954336
63 -5.3003037 -13.1961871
64 7.5852435 -5.3003037
65 7.1698112 7.5852435
66 3.9649379 7.1698112
67 -4.8973156 3.9649379
68 -6.2970058 -4.8973156
69 4.3055179 -6.2970058
70 -8.2650053 4.3055179
71 -22.6019908 -8.2650053
72 -2.2364830 -22.6019908
73 50.0205777 -2.2364830
74 -25.9752327 50.0205777
75 -16.3169407 -25.9752327
76 -6.5873174 -16.3169407
77 21.7915147 -6.5873174
78 58.6052452 21.7915147
79 NA 58.6052452
> dum1 <- dum[2:length(myerror),]
> dum1
lag(myerror, k = 1) myerror
[1,] 68.7263705 -1.8790100
[2,] -9.1766451 68.7263705
[3,] -2.6171986 -9.1766451
[4,] -5.8200672 -2.6171986
[5,] -16.5767531 -5.8200672
[6,] 0.3399743 -16.5767531
[7,] 3.4759765 0.3399743
[8,] -2.0321659 3.4759765
[9,] 13.9321823 -2.0321659
[10,] -7.9112684 13.9321823
[11,] -10.3056632 -7.9112684
[12,] 1.2955652 -10.3056632
[13,] -5.4548847 1.2955652
[14,] -7.3415410 -5.4548847
[15,] -44.9514814 -7.3415410
[16,] -3.4465874 -44.9514814
[17,] -25.2656536 -3.4465874
[18,] -3.7260937 -25.2656536
[19,] -1.0880709 -3.7260937
[20,] -6.3323893 -1.0880709
[21,] 43.4871998 -6.3323893
[22,] 22.4979617 43.4871998
[23,] -4.0750140 22.4979617
[24,] -6.4271977 -4.0750140
[25,] 6.0478784 -6.4271977
[26,] -0.2332449 6.0478784
[27,] 0.1855895 -0.2332449
[28,] -2.2981612 0.1855895
[29,] 1.2908424 -2.2981612
[30,] -10.2247933 1.2908424
[31,] -1.0722031 -10.2247933
[32,] 58.1025208 -1.0722031
[33,] -7.9455529 58.1025208
[34,] -1.8943273 -7.9455529
[35,] 20.9423309 -1.8943273
[36,] 20.3787885 20.9423309
[37,] 25.1978958 20.3787885
[38,] -6.1432142 25.1978958
[39,] -32.4763259 -6.1432142
[40,] 0.9886715 -32.4763259
[41,] -11.7303586 0.9886715
[42,] 0.5878749 -11.7303586
[43,] -9.1390664 0.5878749
[44,] -17.2706721 -9.1390664
[45,] -1.0844117 -17.2706721
[46,] 27.3923124 -1.0844117
[47,] 30.1962195 27.3923124
[48,] -6.6575442 30.1962195
[49,] -5.4868842 -6.6575442
[50,] 9.8292404 -5.4868842
[51,] 19.7954284 9.8292404
[52,] 6.7704649 19.7954284
[53,] 18.4512193 6.7704649
[54,] -8.8951248 18.4512193
[55,] -10.2601667 -8.8951248
[56,] -15.7042332 -10.2601667
[57,] -107.4603053 -15.7042332
[58,] -4.4894928 -107.4603053
[59,] -5.2121801 -4.4894928
[60,] -6.1801919 -5.2121801
[61,] -5.3954336 -6.1801919
[62,] -13.1961871 -5.3954336
[63,] -5.3003037 -13.1961871
[64,] 7.5852435 -5.3003037
[65,] 7.1698112 7.5852435
[66,] 3.9649379 7.1698112
[67,] -4.8973156 3.9649379
[68,] -6.2970058 -4.8973156
[69,] 4.3055179 -6.2970058
[70,] -8.2650053 4.3055179
[71,] -22.6019908 -8.2650053
[72,] -2.2364830 -22.6019908
[73,] 50.0205777 -2.2364830
[74,] -25.9752327 50.0205777
[75,] -16.3169407 -25.9752327
[76,] -6.5873174 -16.3169407
[77,] 21.7915147 -6.5873174
[78,] 58.6052452 21.7915147
> z <- as.data.frame(dum1)
> z
lag(myerror, k = 1) myerror
1 68.7263705 -1.8790100
2 -9.1766451 68.7263705
3 -2.6171986 -9.1766451
4 -5.8200672 -2.6171986
5 -16.5767531 -5.8200672
6 0.3399743 -16.5767531
7 3.4759765 0.3399743
8 -2.0321659 3.4759765
9 13.9321823 -2.0321659
10 -7.9112684 13.9321823
11 -10.3056632 -7.9112684
12 1.2955652 -10.3056632
13 -5.4548847 1.2955652
14 -7.3415410 -5.4548847
15 -44.9514814 -7.3415410
16 -3.4465874 -44.9514814
17 -25.2656536 -3.4465874
18 -3.7260937 -25.2656536
19 -1.0880709 -3.7260937
20 -6.3323893 -1.0880709
21 43.4871998 -6.3323893
22 22.4979617 43.4871998
23 -4.0750140 22.4979617
24 -6.4271977 -4.0750140
25 6.0478784 -6.4271977
26 -0.2332449 6.0478784
27 0.1855895 -0.2332449
28 -2.2981612 0.1855895
29 1.2908424 -2.2981612
30 -10.2247933 1.2908424
31 -1.0722031 -10.2247933
32 58.1025208 -1.0722031
33 -7.9455529 58.1025208
34 -1.8943273 -7.9455529
35 20.9423309 -1.8943273
36 20.3787885 20.9423309
37 25.1978958 20.3787885
38 -6.1432142 25.1978958
39 -32.4763259 -6.1432142
40 0.9886715 -32.4763259
41 -11.7303586 0.9886715
42 0.5878749 -11.7303586
43 -9.1390664 0.5878749
44 -17.2706721 -9.1390664
45 -1.0844117 -17.2706721
46 27.3923124 -1.0844117
47 30.1962195 27.3923124
48 -6.6575442 30.1962195
49 -5.4868842 -6.6575442
50 9.8292404 -5.4868842
51 19.7954284 9.8292404
52 6.7704649 19.7954284
53 18.4512193 6.7704649
54 -8.8951248 18.4512193
55 -10.2601667 -8.8951248
56 -15.7042332 -10.2601667
57 -107.4603053 -15.7042332
58 -4.4894928 -107.4603053
59 -5.2121801 -4.4894928
60 -6.1801919 -5.2121801
61 -5.3954336 -6.1801919
62 -13.1961871 -5.3954336
63 -5.3003037 -13.1961871
64 7.5852435 -5.3003037
65 7.1698112 7.5852435
66 3.9649379 7.1698112
67 -4.8973156 3.9649379
68 -6.2970058 -4.8973156
69 4.3055179 -6.2970058
70 -8.2650053 4.3055179
71 -22.6019908 -8.2650053
72 -2.2364830 -22.6019908
73 50.0205777 -2.2364830
74 -25.9752327 50.0205777
75 -16.3169407 -25.9752327
76 -6.5873174 -16.3169407
77 21.7915147 -6.5873174
78 58.6052452 21.7915147
> 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/738d01355997553.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/8suzj1355997553.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/9l19r1355997553.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/1038fz1355997553.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/11az0o1355997553.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/12959h1355997553.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/13x80v1355997553.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/1413tz1355997553.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/15xphj1355997553.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/167pla1355997553.tab")
+ }
>
> try(system("convert tmp/1f0dj1355997553.ps tmp/1f0dj1355997553.png",intern=TRUE))
character(0)
> try(system("convert tmp/229d61355997553.ps tmp/229d61355997553.png",intern=TRUE))
character(0)
> try(system("convert tmp/3ollp1355997553.ps tmp/3ollp1355997553.png",intern=TRUE))
character(0)
> try(system("convert tmp/405l01355997553.ps tmp/405l01355997553.png",intern=TRUE))
character(0)
> try(system("convert tmp/5z6sw1355997553.ps tmp/5z6sw1355997553.png",intern=TRUE))
character(0)
> try(system("convert tmp/65zc81355997553.ps tmp/65zc81355997553.png",intern=TRUE))
character(0)
> try(system("convert tmp/738d01355997553.ps tmp/738d01355997553.png",intern=TRUE))
character(0)
> try(system("convert tmp/8suzj1355997553.ps tmp/8suzj1355997553.png",intern=TRUE))
character(0)
> try(system("convert tmp/9l19r1355997553.ps tmp/9l19r1355997553.png",intern=TRUE))
character(0)
> try(system("convert tmp/1038fz1355997553.ps tmp/1038fz1355997553.png",intern=TRUE))
character(0)
>
>
> proc.time()
user system elapsed
7.115 1.312 8.428