R version 2.8.0 (2008-10-20)
Copyright (C) 2008 The R Foundation for Statistical Computing
ISBN 3-900051-07-0
R is free software and comes with ABSOLUTELY NO WARRANTY.
You are welcome to redistribute it under certain conditions.
Type 'license()' or 'licence()' for distribution details.
Natural language support but running in an English locale
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(0
+ ,162556
+ ,1081
+ ,213118
+ ,6282154
+ ,0
+ ,29790
+ ,309
+ ,81767
+ ,4321023
+ ,0
+ ,87550
+ ,458
+ ,153198
+ ,4111912
+ ,1
+ ,84738
+ ,588
+ ,-26007
+ ,223193
+ ,0
+ ,54660
+ ,302
+ ,126942
+ ,1491348
+ ,0
+ ,42634
+ ,156
+ ,157214
+ ,1629616
+ ,1
+ ,40949
+ ,481
+ ,129352
+ ,1398893
+ ,0
+ ,45187
+ ,353
+ ,234817
+ ,1926517
+ ,0
+ ,37704
+ ,452
+ ,60448
+ ,983660
+ ,0
+ ,16275
+ ,109
+ ,47818
+ ,1443586
+ ,1
+ ,25830
+ ,115
+ ,245546
+ ,1073089
+ ,1
+ ,12679
+ ,110
+ ,48020
+ ,984885
+ ,0
+ ,18014
+ ,239
+ ,-1710
+ ,1405225
+ ,1
+ ,43556
+ ,247
+ ,32648
+ ,227132
+ ,0
+ ,24811
+ ,505
+ ,95350
+ ,929118
+ ,1
+ ,6575
+ ,159
+ ,151352
+ ,1071292
+ ,1
+ ,7123
+ ,109
+ ,288170
+ ,638830
+ ,0
+ ,21950
+ ,519
+ ,114337
+ ,856956
+ ,0
+ ,37597
+ ,248
+ ,37884
+ ,992426
+ ,1
+ ,17821
+ ,373
+ ,122844
+ ,444477
+ ,0
+ ,12988
+ ,119
+ ,82340
+ ,857217
+ ,0
+ ,22330
+ ,84
+ ,79801
+ ,711969
+ ,1
+ ,13326
+ ,102
+ ,165548
+ ,702380
+ ,1
+ ,16189
+ ,295
+ ,116384
+ ,358589
+ ,1
+ ,7146
+ ,105
+ ,134028
+ ,297978
+ ,1
+ ,15824
+ ,64
+ ,63838
+ ,585715
+ ,0
+ ,27664
+ ,282
+ ,74996
+ ,657954
+ ,1
+ ,11920
+ ,182
+ ,31080
+ ,209458
+ ,1
+ ,8568
+ ,37
+ ,32168
+ ,786690
+ ,1
+ ,14416
+ ,361
+ ,49857
+ ,439798
+ ,0
+ ,3369
+ ,28
+ ,87161
+ ,688779
+ ,0
+ ,11819
+ ,85
+ ,106113
+ ,574339
+ ,0
+ ,6984
+ ,45
+ ,80570
+ ,741409
+ ,0
+ ,4519
+ ,49
+ ,102129
+ ,597793
+ ,1
+ ,2220
+ ,22
+ ,301670
+ ,644190
+ ,1
+ ,18562
+ ,155
+ ,102313
+ ,377934
+ ,1
+ ,10327
+ ,91
+ ,88577
+ ,640273
+ ,0
+ ,5336
+ ,81
+ ,112477
+ ,697458
+ ,0
+ ,2365
+ ,79
+ ,191778
+ ,550608
+ ,1
+ ,4069
+ ,145
+ ,79804
+ ,207393
+ ,1
+ ,8636
+ ,855
+ ,128294
+ ,301607
+ ,1
+ ,13718
+ ,61
+ ,96448
+ ,345783
+ ,1
+ ,4525
+ ,226
+ ,93811
+ ,501749
+ ,1
+ ,6869
+ ,105
+ ,117520
+ ,379983
+ ,1
+ ,4628
+ ,62
+ ,69159
+ ,387475
+ ,0
+ ,3689
+ ,25
+ ,101792
+ ,377305
+ ,0
+ ,4891
+ ,217
+ ,210568
+ ,370837
+ ,0
+ ,7489
+ ,322
+ ,136996
+ ,430866
+ ,1
+ ,4901
+ ,84
+ ,121920
+ ,469107
+ ,1
+ ,2284
+ ,33
+ ,76403
+ ,194493
+ ,0
+ ,3160
+ ,108
+ ,108094
+ ,530670
+ ,0
+ ,4150
+ ,150
+ ,134759
+ ,518365
+ ,0
+ ,7285
+ ,115
+ ,188873
+ ,491303
+ ,0
+ ,1134
+ ,162
+ ,146216
+ ,527021
+ ,0
+ ,4658
+ ,158
+ ,156608
+ ,233773
+ ,1
+ ,2384
+ ,97
+ ,61348
+ ,405972
+ ,1
+ ,3748
+ ,9
+ ,50350
+ ,652925
+ ,1
+ ,5371
+ ,66
+ ,87720
+ ,446211
+ ,1
+ ,1285
+ ,107
+ ,99489
+ ,341340
+ ,0
+ ,9327
+ ,101
+ ,87419
+ ,387699
+ ,0
+ ,5565
+ ,47
+ ,94355
+ ,493408
+ ,1
+ ,1528
+ ,38
+ ,60326
+ ,146494
+ ,0
+ ,3122
+ ,34
+ ,94670
+ ,414462
+ ,0
+ ,7561
+ ,87
+ ,82425
+ ,364304
+ ,1
+ ,2675
+ ,79
+ ,59017
+ ,355178
+ ,1
+ ,13253
+ ,947
+ ,90829
+ ,357760
+ ,1
+ ,880
+ ,74
+ ,80791
+ ,261216
+ ,0
+ ,2053
+ ,53
+ ,100423
+ ,397144
+ ,1
+ ,1424
+ ,94
+ ,131116
+ ,374943
+ ,0
+ ,4036
+ ,63
+ ,100269
+ ,424898
+ ,0
+ ,3045
+ ,58
+ ,27330
+ ,202055
+ ,1
+ ,5119
+ ,49
+ ,39039
+ ,378525
+ ,1
+ ,1431
+ ,34
+ ,106885
+ ,310768
+ ,1
+ ,554
+ ,11
+ ,79285
+ ,325738
+ ,1
+ ,1975
+ ,35
+ ,118881
+ ,394510
+ ,0
+ ,1765
+ ,20
+ ,77623
+ ,247060
+ ,1
+ ,1012
+ ,47
+ ,114768
+ ,368078
+ ,1
+ ,810
+ ,43
+ ,74015
+ ,236761
+ ,1
+ ,1280
+ ,117
+ ,69465
+ ,312378
+ ,0
+ ,666
+ ,171
+ ,117869
+ ,339836
+ ,1
+ ,1380
+ ,26
+ ,60982
+ ,347385
+ ,0
+ ,4677
+ ,75
+ ,90131
+ ,426280
+ ,1
+ ,876
+ ,59
+ ,138971
+ ,352850
+ ,1
+ ,814
+ ,18
+ ,39625
+ ,301881
+ ,1
+ ,514
+ ,15
+ ,102725
+ ,377516
+ ,0
+ ,5692
+ ,72
+ ,64239
+ ,357312
+ ,1
+ ,3642
+ ,86
+ ,90262
+ ,458343
+ ,1
+ ,540
+ ,14
+ ,103960
+ ,354228
+ ,1
+ ,2099
+ ,64
+ ,106611
+ ,308636
+ ,1
+ ,567
+ ,11
+ ,103345
+ ,386212
+ ,1
+ ,2001
+ ,52
+ ,95551
+ ,393343
+ ,0
+ ,2949
+ ,41
+ ,82903
+ ,378509
+ ,1
+ ,2253
+ ,99
+ ,63593
+ ,452469
+ ,0
+ ,6533
+ ,75
+ ,126910
+ ,364839
+ ,1
+ ,1889
+ ,45
+ ,37527
+ ,358649
+ ,0
+ ,3055
+ ,43
+ ,60247
+ ,376641
+ ,1
+ ,272
+ ,8
+ ,112995
+ ,429112
+ ,0
+ ,1414
+ ,198
+ ,70184
+ ,330546
+ ,1
+ ,2564
+ ,22
+ ,130140
+ ,403560
+ ,0
+ ,1383
+ ,11
+ ,73221
+ ,317892)
+ ,dim=c(5
+ ,100)
+ ,dimnames=list(c('Group'
+ ,'Costs'
+ ,'Trades'
+ ,'Dividends'
+ ,'Wealth')
+ ,1:100))
> y <- array(NA,dim=c(5,100),dimnames=list(c('Group','Costs','Trades','Dividends','Wealth'),1:100))
> 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 = '5'
> #'GNU S' R Code compiled by R2WASP v. 1.0.44 ()
> #Author: Prof. Dr. P. Wessa
> #To cite this work: AUTHOR(S), (YEAR), YOUR SOFTWARE TITLE (vNUMBER) in Free Statistics Software (v$_version), Office for Research Development and Education, URL http://www.wessa.net/rwasp_YOURPAGE.wasp/
> #Source of accompanying publication: Office for Research, Development, and Education
> #Technical description: Write here your technical program description (don't use hard returns!)
> library(lattice)
> library(lmtest)
Loading required package: zoo
Attaching package: 'zoo'
The following object(s) are masked from package:base :
as.Date.numeric
> n25 <- 25 #minimum number of obs. for Goldfeld-Quandt test
> par1 <- as.numeric(par1)
> x <- t(y)
> k <- length(x[1,])
> n <- length(x[,1])
> x1 <- cbind(x[,par1], x[,1:k!=par1])
> mycolnames <- c(colnames(x)[par1], colnames(x)[1:k!=par1])
> colnames(x1) <- mycolnames #colnames(x)[par1]
> x <- x1
> if (par3 == 'First Differences'){
+ x2 <- array(0, dim=c(n-1,k), dimnames=list(1:(n-1), paste('(1-B)',colnames(x),sep='')))
+ for (i in 1:n-1) {
+ for (j in 1:k) {
+ x2[i,j] <- x[i+1,j] - x[i,j]
+ }
+ }
+ x <- x2
+ }
> if (par2 == 'Include Monthly Dummies'){
+ x2 <- array(0, dim=c(n,11), dimnames=list(1:n, paste('M', seq(1:11), sep ='')))
+ for (i in 1:11){
+ x2[seq(i,n,12),i] <- 1
+ }
+ x <- cbind(x, x2)
+ }
> if (par2 == 'Include Quarterly Dummies'){
+ x2 <- array(0, dim=c(n,3), dimnames=list(1:n, paste('Q', seq(1:3), sep ='')))
+ for (i in 1:3){
+ x2[seq(i,n,4),i] <- 1
+ }
+ x <- cbind(x, x2)
+ }
> k <- length(x[1,])
> if (par3 == 'Linear Trend'){
+ x <- cbind(x, c(1:n))
+ colnames(x)[k+1] <- 't'
+ }
> x
Wealth Group Costs Trades Dividends
1 6282154 0 162556 1081 213118
2 4321023 0 29790 309 81767
3 4111912 0 87550 458 153198
4 223193 1 84738 588 -26007
5 1491348 0 54660 302 126942
6 1629616 0 42634 156 157214
7 1398893 1 40949 481 129352
8 1926517 0 45187 353 234817
9 983660 0 37704 452 60448
10 1443586 0 16275 109 47818
11 1073089 1 25830 115 245546
12 984885 1 12679 110 48020
13 1405225 0 18014 239 -1710
14 227132 1 43556 247 32648
15 929118 0 24811 505 95350
16 1071292 1 6575 159 151352
17 638830 1 7123 109 288170
18 856956 0 21950 519 114337
19 992426 0 37597 248 37884
20 444477 1 17821 373 122844
21 857217 0 12988 119 82340
22 711969 0 22330 84 79801
23 702380 1 13326 102 165548
24 358589 1 16189 295 116384
25 297978 1 7146 105 134028
26 585715 1 15824 64 63838
27 657954 0 27664 282 74996
28 209458 1 11920 182 31080
29 786690 1 8568 37 32168
30 439798 1 14416 361 49857
31 688779 0 3369 28 87161
32 574339 0 11819 85 106113
33 741409 0 6984 45 80570
34 597793 0 4519 49 102129
35 644190 1 2220 22 301670
36 377934 1 18562 155 102313
37 640273 1 10327 91 88577
38 697458 0 5336 81 112477
39 550608 0 2365 79 191778
40 207393 1 4069 145 79804
41 301607 1 8636 855 128294
42 345783 1 13718 61 96448
43 501749 1 4525 226 93811
44 379983 1 6869 105 117520
45 387475 1 4628 62 69159
46 377305 0 3689 25 101792
47 370837 0 4891 217 210568
48 430866 0 7489 322 136996
49 469107 1 4901 84 121920
50 194493 1 2284 33 76403
51 530670 0 3160 108 108094
52 518365 0 4150 150 134759
53 491303 0 7285 115 188873
54 527021 0 1134 162 146216
55 233773 0 4658 158 156608
56 405972 1 2384 97 61348
57 652925 1 3748 9 50350
58 446211 1 5371 66 87720
59 341340 1 1285 107 99489
60 387699 0 9327 101 87419
61 493408 0 5565 47 94355
62 146494 1 1528 38 60326
63 414462 0 3122 34 94670
64 364304 0 7561 87 82425
65 355178 1 2675 79 59017
66 357760 1 13253 947 90829
67 261216 1 880 74 80791
68 397144 0 2053 53 100423
69 374943 1 1424 94 131116
70 424898 0 4036 63 100269
71 202055 0 3045 58 27330
72 378525 1 5119 49 39039
73 310768 1 1431 34 106885
74 325738 1 554 11 79285
75 394510 1 1975 35 118881
76 247060 0 1765 20 77623
77 368078 1 1012 47 114768
78 236761 1 810 43 74015
79 312378 1 1280 117 69465
80 339836 0 666 171 117869
81 347385 1 1380 26 60982
82 426280 0 4677 75 90131
83 352850 1 876 59 138971
84 301881 1 814 18 39625
85 377516 1 514 15 102725
86 357312 0 5692 72 64239
87 458343 1 3642 86 90262
88 354228 1 540 14 103960
89 308636 1 2099 64 106611
90 386212 1 567 11 103345
91 393343 1 2001 52 95551
92 378509 0 2949 41 82903
93 452469 1 2253 99 63593
94 364839 0 6533 75 126910
95 358649 1 1889 45 37527
96 376641 0 3055 43 60247
97 429112 1 272 8 112995
98 330546 0 1414 198 70184
99 403560 1 2564 22 130140
100 317892 0 1383 11 73221
> k <- length(x[1,])
> df <- as.data.frame(x)
> (mylm <- lm(df))
Call:
lm(formula = df)
Coefficients:
(Intercept) Group Costs Trades Dividends
2.052e+05 -2.103e+05 3.068e+01 -3.246e+02 2.373e+00
> (mysum <- summary(mylm))
Call:
lm(formula = df)
Residuals:
Min 1Q Median 3Q Max
-2118699 -167993 -13843 128253 3108199
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 2.052e+05 1.266e+05 1.620 0.1084
Group -2.103e+05 1.002e+05 -2.098 0.0386 *
Costs 3.068e+01 3.191e+00 9.615 1.1e-15 ***
Trades -3.246e+02 3.586e+02 -0.905 0.3677
Dividends 2.373e+00 9.269e-01 2.560 0.0120 *
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Residual standard error: 487500 on 95 degrees of freedom
Multiple R-squared: 0.674, Adjusted R-squared: 0.6603
F-statistic: 49.11 on 4 and 95 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,] 1.0000000 1.149510e-15 5.747548e-16
[2,] 1.0000000 6.501924e-24 3.250962e-24
[3,] 1.0000000 9.838985e-26 4.919492e-26
[4,] 1.0000000 1.093317e-25 5.466584e-26
[5,] 1.0000000 5.569551e-28 2.784775e-28
[6,] 1.0000000 1.199069e-31 5.995345e-32
[7,] 1.0000000 6.145062e-34 3.072531e-34
[8,] 1.0000000 6.804677e-35 3.402339e-35
[9,] 1.0000000 1.813588e-38 9.067940e-39
[10,] 1.0000000 1.006113e-37 5.030564e-38
[11,] 1.0000000 4.146459e-38 2.073229e-38
[12,] 1.0000000 6.457182e-38 3.228591e-38
[13,] 1.0000000 4.805932e-37 2.402966e-37
[14,] 1.0000000 1.287432e-37 6.437162e-38
[15,] 1.0000000 3.965084e-37 1.982542e-37
[16,] 1.0000000 1.390437e-36 6.952183e-37
[17,] 1.0000000 6.303613e-36 3.151807e-36
[18,] 1.0000000 2.304466e-35 1.152233e-35
[19,] 1.0000000 1.072559e-34 5.362794e-35
[20,] 1.0000000 2.822984e-34 1.411492e-34
[21,] 1.0000000 4.275950e-34 2.137975e-34
[22,] 1.0000000 5.710565e-36 2.855282e-36
[23,] 1.0000000 3.735881e-35 1.867940e-35
[24,] 1.0000000 1.524535e-35 7.622677e-36
[25,] 1.0000000 6.170279e-35 3.085139e-35
[26,] 1.0000000 6.030711e-36 3.015355e-36
[27,] 1.0000000 8.718315e-36 4.359157e-36
[28,] 1.0000000 4.625266e-35 2.312633e-35
[29,] 1.0000000 1.460461e-34 7.302307e-35
[30,] 1.0000000 1.126734e-34 5.633670e-35
[31,] 1.0000000 8.826448e-36 4.413224e-36
[32,] 1.0000000 2.496123e-35 1.248061e-35
[33,] 1.0000000 4.685161e-35 2.342581e-35
[34,] 1.0000000 3.159456e-34 1.579728e-34
[35,] 1.0000000 8.850116e-34 4.425058e-34
[36,] 1.0000000 2.131976e-33 1.065988e-33
[37,] 1.0000000 1.455305e-32 7.276523e-33
[38,] 1.0000000 1.084196e-31 5.420979e-32
[39,] 1.0000000 7.412173e-31 3.706087e-31
[40,] 1.0000000 2.027283e-30 1.013641e-30
[41,] 1.0000000 1.330207e-29 6.651035e-30
[42,] 1.0000000 8.833899e-29 4.416949e-29
[43,] 1.0000000 1.098503e-28 5.492516e-29
[44,] 1.0000000 1.660361e-28 8.301805e-29
[45,] 1.0000000 4.020852e-28 2.010426e-28
[46,] 1.0000000 2.465482e-27 1.232741e-27
[47,] 1.0000000 8.748733e-28 4.374366e-28
[48,] 1.0000000 8.178200e-28 4.089100e-28
[49,] 1.0000000 4.680325e-27 2.340163e-27
[50,] 1.0000000 5.301790e-30 2.650895e-30
[51,] 1.0000000 4.019791e-29 2.009895e-29
[52,] 1.0000000 4.113444e-28 2.056722e-28
[53,] 1.0000000 2.737689e-27 1.368844e-27
[54,] 1.0000000 9.166718e-27 4.583359e-27
[55,] 1.0000000 5.886556e-28 2.943278e-28
[56,] 1.0000000 4.196695e-27 2.098347e-27
[57,] 1.0000000 3.585171e-26 1.792585e-26
[58,] 1.0000000 3.887099e-25 1.943549e-25
[59,] 1.0000000 2.634617e-24 1.317308e-24
[60,] 1.0000000 9.038535e-24 4.519267e-24
[61,] 1.0000000 5.550277e-23 2.775138e-23
[62,] 1.0000000 5.822142e-22 2.911071e-22
[63,] 1.0000000 3.457048e-21 1.728524e-21
[64,] 1.0000000 4.960305e-21 2.480152e-21
[65,] 1.0000000 5.302941e-20 2.651471e-20
[66,] 1.0000000 3.137319e-19 1.568660e-19
[67,] 1.0000000 3.363604e-18 1.681802e-18
[68,] 1.0000000 3.714660e-17 1.857330e-17
[69,] 1.0000000 1.167859e-16 5.839297e-17
[70,] 1.0000000 1.312036e-15 6.560181e-16
[71,] 1.0000000 6.973451e-16 3.486726e-16
[72,] 1.0000000 3.744969e-15 1.872485e-15
[73,] 1.0000000 4.758931e-14 2.379466e-14
[74,] 1.0000000 5.093667e-13 2.546834e-13
[75,] 1.0000000 3.142222e-12 1.571111e-12
[76,] 1.0000000 2.826536e-11 1.413268e-11
[77,] 1.0000000 9.543906e-11 4.771953e-11
[78,] 1.0000000 1.208650e-09 6.043251e-10
[79,] 1.0000000 1.319321e-08 6.596604e-09
[80,] 1.0000000 8.340256e-08 4.170128e-08
[81,] 0.9999996 8.184116e-07 4.092058e-07
[82,] 0.9999995 1.040498e-06 5.202489e-07
[83,] 0.9999928 1.435868e-05 7.179340e-06
[84,] 0.9999081 1.838548e-04 9.192742e-05
[85,] 0.9990223 1.955347e-03 9.776735e-04
> postscript(file="/var/www/html/freestat/rcomp/tmp/1y4381291392198.ps",horizontal=F,onefile=F,pagecentre=F,paper="special",width=8.3333333333333,height=5.5555555555556)
> plot(x[,1], type='l', main='Actuals and Interpolation', ylab='value of Actuals and Interpolation (dots)', xlab='time or index')
> points(x[,1]-mysum$resid)
> grid()
> dev.off()
null device
1
> postscript(file="/var/www/html/freestat/rcomp/tmp/2rwks1291392198.ps",horizontal=F,onefile=F,pagecentre=F,paper="special",width=8.3333333333333,height=5.5555555555556)
> plot(mysum$resid, type='b', pch=19, main='Residuals', ylab='value of Residuals', xlab='time or index')
> grid()
> dev.off()
null device
1
> postscript(file="/var/www/html/freestat/rcomp/tmp/3rwks1291392198.ps",horizontal=F,onefile=F,pagecentre=F,paper="special",width=8.3333333333333,height=5.5555555555556)
> hist(mysum$resid, main='Residual Histogram', xlab='values of Residuals')
> grid()
> dev.off()
null device
1
> postscript(file="/var/www/html/freestat/rcomp/tmp/4rwks1291392198.ps",horizontal=F,onefile=F,pagecentre=F,paper="special",width=8.3333333333333,height=5.5555555555556)
> densityplot(~mysum$resid,col='black',main='Residual Density Plot', xlab='values of Residuals')
> dev.off()
null device
1
> postscript(file="/var/www/html/freestat/rcomp/tmp/5252d1291392198.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 = 100
Frequency = 1
1 2 3 4 5 6
935301.699 3108199.099 1006023.279 -2118698.544 -593892.459 -205923.000
7 8 9 10 11 12
-3072.106 -107555.022 -374939.913 661013.370 -259589.926 522753.913
13 14 15 16 17 18
729029.653 -1101276.542 -99581.379 567109.219 -223064.876 -124488.452
19 20 21 22 23 24
-375558.353 -267594.388 96804.587 -340365.032 -61092.721 -313403.670
25 26 27 28 29 30
-200137.866 -25362.669 -482345.650 -165821.120 464596.823 1487.822
31 32 33 34 35 36
182475.008 -217660.947 145365.018 27506.656 -127567.573 -378904.444
37 38 39 40 41 42
147885.227 87938.450 -156601.856 -54671.337 14814.867 -279046.440
43 44 45 46 47 48
218746.871 -70461.368 106583.801 -174509.148 -413662.363 -224665.013
49 50 51 52 53 54
61778.221 -41093.392 7068.445 -85251.898 -348260.836 -7366.253
55 56 57 58 59 60
-434679.626 223816.123 426461.914 99779.156 105634.541 -278298.824
61 62 63 64 65 66
-91167.984 -26126.426 -100136.351 -240210.865 163784.343 48096.647
67 68 69 70 71 72
71594.913 -92145.597 55702.352 -121613.456 -142592.054 149827.410
73 74 75 76 77 78
29339.313 129243.696 68250.647 -190000.314 85015.905 55305.336
79 80 81 82 83 84
151319.316 -110005.043 173853.288 -111943.028 20420.495 193796.768
85 86 87 88 89 90
127923.346 -151579.769 165411.246 100582.483 17102.157 132223.904
91 92 93 94 95 96
127166.372 -100586.824 269653.536 -317598.497 231328.665 -51294.217
97 98 99 100
160300.295 -20320.583 28294.474 -99924.659
> postscript(file="/var/www/html/freestat/rcomp/tmp/6252d1291392198.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 = 100
Frequency = 1
lag(myerror, k = 1) myerror
0 935301.699 NA
1 3108199.099 935301.699
2 1006023.279 3108199.099
3 -2118698.544 1006023.279
4 -593892.459 -2118698.544
5 -205923.000 -593892.459
6 -3072.106 -205923.000
7 -107555.022 -3072.106
8 -374939.913 -107555.022
9 661013.370 -374939.913
10 -259589.926 661013.370
11 522753.913 -259589.926
12 729029.653 522753.913
13 -1101276.542 729029.653
14 -99581.379 -1101276.542
15 567109.219 -99581.379
16 -223064.876 567109.219
17 -124488.452 -223064.876
18 -375558.353 -124488.452
19 -267594.388 -375558.353
20 96804.587 -267594.388
21 -340365.032 96804.587
22 -61092.721 -340365.032
23 -313403.670 -61092.721
24 -200137.866 -313403.670
25 -25362.669 -200137.866
26 -482345.650 -25362.669
27 -165821.120 -482345.650
28 464596.823 -165821.120
29 1487.822 464596.823
30 182475.008 1487.822
31 -217660.947 182475.008
32 145365.018 -217660.947
33 27506.656 145365.018
34 -127567.573 27506.656
35 -378904.444 -127567.573
36 147885.227 -378904.444
37 87938.450 147885.227
38 -156601.856 87938.450
39 -54671.337 -156601.856
40 14814.867 -54671.337
41 -279046.440 14814.867
42 218746.871 -279046.440
43 -70461.368 218746.871
44 106583.801 -70461.368
45 -174509.148 106583.801
46 -413662.363 -174509.148
47 -224665.013 -413662.363
48 61778.221 -224665.013
49 -41093.392 61778.221
50 7068.445 -41093.392
51 -85251.898 7068.445
52 -348260.836 -85251.898
53 -7366.253 -348260.836
54 -434679.626 -7366.253
55 223816.123 -434679.626
56 426461.914 223816.123
57 99779.156 426461.914
58 105634.541 99779.156
59 -278298.824 105634.541
60 -91167.984 -278298.824
61 -26126.426 -91167.984
62 -100136.351 -26126.426
63 -240210.865 -100136.351
64 163784.343 -240210.865
65 48096.647 163784.343
66 71594.913 48096.647
67 -92145.597 71594.913
68 55702.352 -92145.597
69 -121613.456 55702.352
70 -142592.054 -121613.456
71 149827.410 -142592.054
72 29339.313 149827.410
73 129243.696 29339.313
74 68250.647 129243.696
75 -190000.314 68250.647
76 85015.905 -190000.314
77 55305.336 85015.905
78 151319.316 55305.336
79 -110005.043 151319.316
80 173853.288 -110005.043
81 -111943.028 173853.288
82 20420.495 -111943.028
83 193796.768 20420.495
84 127923.346 193796.768
85 -151579.769 127923.346
86 165411.246 -151579.769
87 100582.483 165411.246
88 17102.157 100582.483
89 132223.904 17102.157
90 127166.372 132223.904
91 -100586.824 127166.372
92 269653.536 -100586.824
93 -317598.497 269653.536
94 231328.665 -317598.497
95 -51294.217 231328.665
96 160300.295 -51294.217
97 -20320.583 160300.295
98 28294.474 -20320.583
99 -99924.659 28294.474
100 NA -99924.659
> dum1 <- dum[2:length(myerror),]
> dum1
lag(myerror, k = 1) myerror
[1,] 3108199.099 935301.699
[2,] 1006023.279 3108199.099
[3,] -2118698.544 1006023.279
[4,] -593892.459 -2118698.544
[5,] -205923.000 -593892.459
[6,] -3072.106 -205923.000
[7,] -107555.022 -3072.106
[8,] -374939.913 -107555.022
[9,] 661013.370 -374939.913
[10,] -259589.926 661013.370
[11,] 522753.913 -259589.926
[12,] 729029.653 522753.913
[13,] -1101276.542 729029.653
[14,] -99581.379 -1101276.542
[15,] 567109.219 -99581.379
[16,] -223064.876 567109.219
[17,] -124488.452 -223064.876
[18,] -375558.353 -124488.452
[19,] -267594.388 -375558.353
[20,] 96804.587 -267594.388
[21,] -340365.032 96804.587
[22,] -61092.721 -340365.032
[23,] -313403.670 -61092.721
[24,] -200137.866 -313403.670
[25,] -25362.669 -200137.866
[26,] -482345.650 -25362.669
[27,] -165821.120 -482345.650
[28,] 464596.823 -165821.120
[29,] 1487.822 464596.823
[30,] 182475.008 1487.822
[31,] -217660.947 182475.008
[32,] 145365.018 -217660.947
[33,] 27506.656 145365.018
[34,] -127567.573 27506.656
[35,] -378904.444 -127567.573
[36,] 147885.227 -378904.444
[37,] 87938.450 147885.227
[38,] -156601.856 87938.450
[39,] -54671.337 -156601.856
[40,] 14814.867 -54671.337
[41,] -279046.440 14814.867
[42,] 218746.871 -279046.440
[43,] -70461.368 218746.871
[44,] 106583.801 -70461.368
[45,] -174509.148 106583.801
[46,] -413662.363 -174509.148
[47,] -224665.013 -413662.363
[48,] 61778.221 -224665.013
[49,] -41093.392 61778.221
[50,] 7068.445 -41093.392
[51,] -85251.898 7068.445
[52,] -348260.836 -85251.898
[53,] -7366.253 -348260.836
[54,] -434679.626 -7366.253
[55,] 223816.123 -434679.626
[56,] 426461.914 223816.123
[57,] 99779.156 426461.914
[58,] 105634.541 99779.156
[59,] -278298.824 105634.541
[60,] -91167.984 -278298.824
[61,] -26126.426 -91167.984
[62,] -100136.351 -26126.426
[63,] -240210.865 -100136.351
[64,] 163784.343 -240210.865
[65,] 48096.647 163784.343
[66,] 71594.913 48096.647
[67,] -92145.597 71594.913
[68,] 55702.352 -92145.597
[69,] -121613.456 55702.352
[70,] -142592.054 -121613.456
[71,] 149827.410 -142592.054
[72,] 29339.313 149827.410
[73,] 129243.696 29339.313
[74,] 68250.647 129243.696
[75,] -190000.314 68250.647
[76,] 85015.905 -190000.314
[77,] 55305.336 85015.905
[78,] 151319.316 55305.336
[79,] -110005.043 151319.316
[80,] 173853.288 -110005.043
[81,] -111943.028 173853.288
[82,] 20420.495 -111943.028
[83,] 193796.768 20420.495
[84,] 127923.346 193796.768
[85,] -151579.769 127923.346
[86,] 165411.246 -151579.769
[87,] 100582.483 165411.246
[88,] 17102.157 100582.483
[89,] 132223.904 17102.157
[90,] 127166.372 132223.904
[91,] -100586.824 127166.372
[92,] 269653.536 -100586.824
[93,] -317598.497 269653.536
[94,] 231328.665 -317598.497
[95,] -51294.217 231328.665
[96,] 160300.295 -51294.217
[97,] -20320.583 160300.295
[98,] 28294.474 -20320.583
[99,] -99924.659 28294.474
> z <- as.data.frame(dum1)
> z
lag(myerror, k = 1) myerror
1 3108199.099 935301.699
2 1006023.279 3108199.099
3 -2118698.544 1006023.279
4 -593892.459 -2118698.544
5 -205923.000 -593892.459
6 -3072.106 -205923.000
7 -107555.022 -3072.106
8 -374939.913 -107555.022
9 661013.370 -374939.913
10 -259589.926 661013.370
11 522753.913 -259589.926
12 729029.653 522753.913
13 -1101276.542 729029.653
14 -99581.379 -1101276.542
15 567109.219 -99581.379
16 -223064.876 567109.219
17 -124488.452 -223064.876
18 -375558.353 -124488.452
19 -267594.388 -375558.353
20 96804.587 -267594.388
21 -340365.032 96804.587
22 -61092.721 -340365.032
23 -313403.670 -61092.721
24 -200137.866 -313403.670
25 -25362.669 -200137.866
26 -482345.650 -25362.669
27 -165821.120 -482345.650
28 464596.823 -165821.120
29 1487.822 464596.823
30 182475.008 1487.822
31 -217660.947 182475.008
32 145365.018 -217660.947
33 27506.656 145365.018
34 -127567.573 27506.656
35 -378904.444 -127567.573
36 147885.227 -378904.444
37 87938.450 147885.227
38 -156601.856 87938.450
39 -54671.337 -156601.856
40 14814.867 -54671.337
41 -279046.440 14814.867
42 218746.871 -279046.440
43 -70461.368 218746.871
44 106583.801 -70461.368
45 -174509.148 106583.801
46 -413662.363 -174509.148
47 -224665.013 -413662.363
48 61778.221 -224665.013
49 -41093.392 61778.221
50 7068.445 -41093.392
51 -85251.898 7068.445
52 -348260.836 -85251.898
53 -7366.253 -348260.836
54 -434679.626 -7366.253
55 223816.123 -434679.626
56 426461.914 223816.123
57 99779.156 426461.914
58 105634.541 99779.156
59 -278298.824 105634.541
60 -91167.984 -278298.824
61 -26126.426 -91167.984
62 -100136.351 -26126.426
63 -240210.865 -100136.351
64 163784.343 -240210.865
65 48096.647 163784.343
66 71594.913 48096.647
67 -92145.597 71594.913
68 55702.352 -92145.597
69 -121613.456 55702.352
70 -142592.054 -121613.456
71 149827.410 -142592.054
72 29339.313 149827.410
73 129243.696 29339.313
74 68250.647 129243.696
75 -190000.314 68250.647
76 85015.905 -190000.314
77 55305.336 85015.905
78 151319.316 55305.336
79 -110005.043 151319.316
80 173853.288 -110005.043
81 -111943.028 173853.288
82 20420.495 -111943.028
83 193796.768 20420.495
84 127923.346 193796.768
85 -151579.769 127923.346
86 165411.246 -151579.769
87 100582.483 165411.246
88 17102.157 100582.483
89 132223.904 17102.157
90 127166.372 132223.904
91 -100586.824 127166.372
92 269653.536 -100586.824
93 -317598.497 269653.536
94 231328.665 -317598.497
95 -51294.217 231328.665
96 160300.295 -51294.217
97 -20320.583 160300.295
98 28294.474 -20320.583
99 -99924.659 28294.474
> plot(z,main=paste('Residual Lag plot, lowess, and regression line'), ylab='values of Residuals', xlab='lagged values of Residuals')
> lines(lowess(z))
> abline(lm(z))
> grid()
> dev.off()
null device
1
> postscript(file="/var/www/html/freestat/rcomp/tmp/7cejg1291392198.ps",horizontal=F,onefile=F,pagecentre=F,paper="special",width=8.3333333333333,height=5.5555555555556)
> acf(mysum$resid, lag.max=length(mysum$resid)/2, main='Residual Autocorrelation Function')
> grid()
> dev.off()
null device
1
> postscript(file="/var/www/html/freestat/rcomp/tmp/8cejg1291392198.ps",horizontal=F,onefile=F,pagecentre=F,paper="special",width=8.3333333333333,height=5.5555555555556)
> pacf(mysum$resid, lag.max=length(mysum$resid)/2, main='Residual Partial Autocorrelation Function')
> grid()
> dev.off()
null device
1
> postscript(file="/var/www/html/freestat/rcomp/tmp/9nn011291392198.ps",horizontal=F,onefile=F,pagecentre=F,paper="special",width=8.3333333333333,height=5.5555555555556)
> opar <- par(mfrow = c(2,2), oma = c(0, 0, 1.1, 0))
> plot(mylm, las = 1, sub='Residual Diagnostics')
> par(opar)
> dev.off()
null device
1
> if (n > n25) {
+ postscript(file="/var/www/html/freestat/rcomp/tmp/10nn011291392198.ps",horizontal=F,onefile=F,pagecentre=F,paper="special",width=8.3333333333333,height=5.5555555555556)
+ plot(kp3:nmkm3,gqarr[,2], main='Goldfeld-Quandt test',ylab='2-sided p-value',xlab='breakpoint')
+ grid()
+ dev.off()
+ }
null device
1
>
> #Note: the /var/www/html/freestat/rcomp/createtable file can be downloaded at http://www.wessa.net/cretab
> load(file="/var/www/html/freestat/rcomp/createtable")
>
> a<-table.start()
> a<-table.row.start(a)
> a<-table.element(a, 'Multiple Linear Regression - Estimated Regression Equation', 1, TRUE)
> a<-table.row.end(a)
> myeq <- colnames(x)[1]
> myeq <- paste(myeq, '[t] = ', sep='')
> for (i in 1:k){
+ if (mysum$coefficients[i,1] > 0) myeq <- paste(myeq, '+', '')
+ myeq <- paste(myeq, mysum$coefficients[i,1], sep=' ')
+ if (rownames(mysum$coefficients)[i] != '(Intercept)') {
+ myeq <- paste(myeq, rownames(mysum$coefficients)[i], sep='')
+ if (rownames(mysum$coefficients)[i] != 't') myeq <- paste(myeq, '[t]', sep='')
+ }
+ }
> myeq <- paste(myeq, ' + e[t]')
> a<-table.row.start(a)
> a<-table.element(a, myeq)
> a<-table.row.end(a)
> a<-table.end(a)
> table.save(a,file="/var/www/html/freestat/rcomp/tmp/1186hp1291392198.tab")
> a<-table.start()
> a<-table.row.start(a)
> a<-table.element(a,hyperlink('http://www.xycoon.com/ols1.htm','Multiple Linear Regression - Ordinary Least Squares',''), 6, TRUE)
> a<-table.row.end(a)
> a<-table.row.start(a)
> a<-table.element(a,'Variable',header=TRUE)
> a<-table.element(a,'Parameter',header=TRUE)
> a<-table.element(a,'S.D.',header=TRUE)
> a<-table.element(a,'T-STAT
H0: parameter = 0',header=TRUE)
> a<-table.element(a,'2-tail p-value',header=TRUE)
> a<-table.element(a,'1-tail p-value',header=TRUE)
> a<-table.row.end(a)
> for (i in 1:k){
+ a<-table.row.start(a)
+ a<-table.element(a,rownames(mysum$coefficients)[i],header=TRUE)
+ a<-table.element(a,mysum$coefficients[i,1])
+ a<-table.element(a, round(mysum$coefficients[i,2],6))
+ a<-table.element(a, round(mysum$coefficients[i,3],4))
+ a<-table.element(a, round(mysum$coefficients[i,4],6))
+ a<-table.element(a, round(mysum$coefficients[i,4]/2,6))
+ a<-table.row.end(a)
+ }
> a<-table.end(a)
> table.save(a,file="/var/www/html/freestat/rcomp/tmp/12coxd1291392198.tab")
> a<-table.start()
> a<-table.row.start(a)
> a<-table.element(a, 'Multiple Linear Regression - Regression Statistics', 2, TRUE)
> a<-table.row.end(a)
> a<-table.row.start(a)
> a<-table.element(a, 'Multiple R',1,TRUE)
> a<-table.element(a, sqrt(mysum$r.squared))
> a<-table.row.end(a)
> a<-table.row.start(a)
> a<-table.element(a, 'R-squared',1,TRUE)
> a<-table.element(a, mysum$r.squared)
> a<-table.row.end(a)
> a<-table.row.start(a)
> a<-table.element(a, 'Adjusted R-squared',1,TRUE)
> a<-table.element(a, mysum$adj.r.squared)
> a<-table.row.end(a)
> a<-table.row.start(a)
> a<-table.element(a, 'F-TEST (value)',1,TRUE)
> a<-table.element(a, mysum$fstatistic[1])
> a<-table.row.end(a)
> a<-table.row.start(a)
> a<-table.element(a, 'F-TEST (DF numerator)',1,TRUE)
> a<-table.element(a, mysum$fstatistic[2])
> a<-table.row.end(a)
> a<-table.row.start(a)
> a<-table.element(a, 'F-TEST (DF denominator)',1,TRUE)
> a<-table.element(a, mysum$fstatistic[3])
> a<-table.row.end(a)
> a<-table.row.start(a)
> a<-table.element(a, 'p-value',1,TRUE)
> a<-table.element(a, 1-pf(mysum$fstatistic[1],mysum$fstatistic[2],mysum$fstatistic[3]))
> a<-table.row.end(a)
> a<-table.row.start(a)
> a<-table.element(a, 'Multiple Linear Regression - Residual Statistics', 2, TRUE)
> a<-table.row.end(a)
> a<-table.row.start(a)
> a<-table.element(a, 'Residual Standard Deviation',1,TRUE)
> a<-table.element(a, mysum$sigma)
> a<-table.row.end(a)
> a<-table.row.start(a)
> a<-table.element(a, 'Sum Squared Residuals',1,TRUE)
> a<-table.element(a, sum(myerror*myerror))
> a<-table.row.end(a)
> a<-table.end(a)
> table.save(a,file="/var/www/html/freestat/rcomp/tmp/131qcp1291392198.tab")
> a<-table.start()
> a<-table.row.start(a)
> a<-table.element(a, 'Multiple Linear Regression - Actuals, Interpolation, and Residuals', 4, TRUE)
> a<-table.row.end(a)
> a<-table.row.start(a)
> a<-table.element(a, 'Time or Index', 1, TRUE)
> a<-table.element(a, 'Actuals', 1, TRUE)
> a<-table.element(a, 'Interpolation
Forecast', 1, TRUE)
> a<-table.element(a, 'Residuals
Prediction Error', 1, TRUE)
> a<-table.row.end(a)
> for (i in 1:n) {
+ a<-table.row.start(a)
+ a<-table.element(a,i, 1, TRUE)
+ a<-table.element(a,x[i])
+ a<-table.element(a,x[i]-mysum$resid[i])
+ a<-table.element(a,mysum$resid[i])
+ a<-table.row.end(a)
+ }
> a<-table.end(a)
> table.save(a,file="/var/www/html/freestat/rcomp/tmp/14thtr1291392198.tab")
> if (n > n25) {
+ a<-table.start()
+ a<-table.row.start(a)
+ a<-table.element(a,'Goldfeld-Quandt test for Heteroskedasticity',4,TRUE)
+ a<-table.row.end(a)
+ a<-table.row.start(a)
+ a<-table.element(a,'p-values',header=TRUE)
+ a<-table.element(a,'Alternative Hypothesis',3,header=TRUE)
+ a<-table.row.end(a)
+ a<-table.row.start(a)
+ a<-table.element(a,'breakpoint index',header=TRUE)
+ a<-table.element(a,'greater',header=TRUE)
+ a<-table.element(a,'2-sided',header=TRUE)
+ a<-table.element(a,'less',header=TRUE)
+ a<-table.row.end(a)
+ for (mypoint in kp3:nmkm3) {
+ a<-table.row.start(a)
+ a<-table.element(a,mypoint,header=TRUE)
+ a<-table.element(a,gqarr[mypoint-kp3+1,1])
+ a<-table.element(a,gqarr[mypoint-kp3+1,2])
+ a<-table.element(a,gqarr[mypoint-kp3+1,3])
+ a<-table.row.end(a)
+ }
+ a<-table.end(a)
+ table.save(a,file="/var/www/html/freestat/rcomp/tmp/15ezaf1291392198.tab")
+ a<-table.start()
+ a<-table.row.start(a)
+ a<-table.element(a,'Meta Analysis of Goldfeld-Quandt test for Heteroskedasticity',4,TRUE)
+ a<-table.row.end(a)
+ a<-table.row.start(a)
+ a<-table.element(a,'Description',header=TRUE)
+ a<-table.element(a,'# significant tests',header=TRUE)
+ a<-table.element(a,'% significant tests',header=TRUE)
+ a<-table.element(a,'OK/NOK',header=TRUE)
+ a<-table.row.end(a)
+ a<-table.row.start(a)
+ a<-table.element(a,'1% type I error level',header=TRUE)
+ a<-table.element(a,numsignificant1)
+ a<-table.element(a,numsignificant1/numgqtests)
+ if (numsignificant1/numgqtests < 0.01) dum <- 'OK' else dum <- 'NOK'
+ a<-table.element(a,dum)
+ a<-table.row.end(a)
+ a<-table.row.start(a)
+ a<-table.element(a,'5% type I error level',header=TRUE)
+ a<-table.element(a,numsignificant5)
+ a<-table.element(a,numsignificant5/numgqtests)
+ if (numsignificant5/numgqtests < 0.05) dum <- 'OK' else dum <- 'NOK'
+ a<-table.element(a,dum)
+ a<-table.row.end(a)
+ a<-table.row.start(a)
+ a<-table.element(a,'10% type I error level',header=TRUE)
+ a<-table.element(a,numsignificant10)
+ a<-table.element(a,numsignificant10/numgqtests)
+ if (numsignificant10/numgqtests < 0.1) dum <- 'OK' else dum <- 'NOK'
+ a<-table.element(a,dum)
+ a<-table.row.end(a)
+ a<-table.end(a)
+ table.save(a,file="/var/www/html/freestat/rcomp/tmp/16t98o1291392198.tab")
+ }
>
> try(system("convert tmp/1y4381291392198.ps tmp/1y4381291392198.png",intern=TRUE))
character(0)
> try(system("convert tmp/2rwks1291392198.ps tmp/2rwks1291392198.png",intern=TRUE))
character(0)
> try(system("convert tmp/3rwks1291392198.ps tmp/3rwks1291392198.png",intern=TRUE))
character(0)
> try(system("convert tmp/4rwks1291392198.ps tmp/4rwks1291392198.png",intern=TRUE))
character(0)
> try(system("convert tmp/5252d1291392198.ps tmp/5252d1291392198.png",intern=TRUE))
character(0)
> try(system("convert tmp/6252d1291392198.ps tmp/6252d1291392198.png",intern=TRUE))
character(0)
> try(system("convert tmp/7cejg1291392198.ps tmp/7cejg1291392198.png",intern=TRUE))
character(0)
> try(system("convert tmp/8cejg1291392198.ps tmp/8cejg1291392198.png",intern=TRUE))
character(0)
> try(system("convert tmp/9nn011291392198.ps tmp/9nn011291392198.png",intern=TRUE))
character(0)
> try(system("convert tmp/10nn011291392198.ps tmp/10nn011291392198.png",intern=TRUE))
character(0)
>
>
> proc.time()
user system elapsed
4.589 2.671 5.066