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(1
+ ,21221
+ ,6086
+ ,7
+ ,2
+ ,7106
+ ,2273
+ ,6
+ ,3
+ ,4794
+ ,2601
+ ,6
+ ,4
+ ,6776
+ ,2708
+ ,6
+ ,5
+ ,7122
+ ,3321
+ ,6
+ ,6
+ ,12025
+ ,2106
+ ,6
+ ,7
+ ,18945
+ ,6353
+ ,5
+ ,8
+ ,24380
+ ,7912
+ ,5
+ ,9
+ ,6942
+ ,1905
+ ,5
+ ,10
+ ,7370
+ ,3841
+ ,5
+ ,11
+ ,13143
+ ,2737
+ ,5
+ ,12
+ ,11831
+ ,4449
+ ,5
+ ,13
+ ,3102
+ ,1253
+ ,5
+ ,14
+ ,7355
+ ,2372
+ ,5
+ ,15
+ ,7749
+ ,3303
+ ,5
+ ,16
+ ,10306
+ ,2876
+ ,5
+ ,17
+ ,7860
+ ,2939
+ ,4
+ ,18
+ ,9138
+ ,2149
+ ,4
+ ,19
+ ,6545
+ ,2545
+ ,4
+ ,20
+ ,9829
+ ,3734
+ ,4
+ ,21
+ ,9460
+ ,1911
+ ,4
+ ,22
+ ,30631
+ ,875
+ ,4
+ ,23
+ ,4807
+ ,1466
+ ,4
+ ,24
+ ,4618
+ ,1562
+ ,4
+ ,25
+ ,7304
+ ,1203
+ ,4
+ ,26
+ ,12621
+ ,2860
+ ,4
+ ,27
+ ,9320
+ ,3245
+ ,4
+ ,28
+ ,9569
+ ,2668
+ ,4
+ ,29
+ ,8480
+ ,2111
+ ,4
+ ,30
+ ,4718
+ ,942
+ ,4
+ ,31
+ ,7342
+ ,2690
+ ,4
+ ,32
+ ,8230
+ ,3675
+ ,4
+ ,33
+ ,10007
+ ,3798
+ ,4
+ ,34
+ ,9176
+ ,672
+ ,4
+ ,35
+ ,2334
+ ,523
+ ,4
+ ,36
+ ,5459
+ ,2352
+ ,3
+ ,37
+ ,5440
+ ,1476
+ ,3
+ ,38
+ ,9027
+ ,2646
+ ,3
+ ,39
+ ,9852
+ ,2500
+ ,3
+ ,40
+ ,3122
+ ,720
+ ,3
+ ,41
+ ,8863
+ ,2483
+ ,3
+ ,42
+ ,4286
+ ,1550
+ ,3
+ ,43
+ ,9119
+ ,2551
+ ,3
+ ,44
+ ,7297
+ ,2544
+ ,3
+ ,45
+ ,8009
+ ,2870
+ ,3
+ ,46
+ ,5049
+ ,1288
+ ,3
+ ,47
+ ,5057
+ ,1474
+ ,3
+ ,48
+ ,4372
+ ,1521
+ ,3
+ ,49
+ ,3393
+ ,1331
+ ,3
+ ,50
+ ,7586
+ ,1685
+ ,3
+ ,51
+ ,7440
+ ,2755
+ ,3
+ ,52
+ ,8636
+ ,1268
+ ,3
+ ,53
+ ,10846
+ ,1266
+ ,3
+ ,54
+ ,6843
+ ,1564
+ ,3
+ ,55
+ ,5102
+ ,1062
+ ,3
+ ,56
+ ,4183
+ ,1306
+ ,3
+ ,57
+ ,10083
+ ,4566
+ ,3
+ ,58
+ ,5273
+ ,1655
+ ,3
+ ,59
+ ,6433
+ ,1423
+ ,3
+ ,60
+ ,6116
+ ,474
+ ,3
+ ,61
+ ,5887
+ ,725
+ ,3
+ ,62
+ ,7415
+ ,1061
+ ,2
+ ,63
+ ,7953
+ ,687
+ ,2
+ ,64
+ ,8447
+ ,2411
+ ,2
+ ,65
+ ,4733
+ ,1705
+ ,2
+ ,66
+ ,11478
+ ,2422
+ ,2
+ ,67
+ ,4234
+ ,1900
+ ,2
+ ,68
+ ,11089
+ ,2848
+ ,2
+ ,69
+ ,6390
+ ,1398
+ ,2
+ ,70
+ ,7737
+ ,2625
+ ,2)
+ ,dim=c(4
+ ,70)
+ ,dimnames=list(c('Ranking'
+ ,'Characters'
+ ,'Revisions'
+ ,'Hours
')
+ ,1:70))
> y <- array(NA,dim=c(4,70),dimnames=list(c('Ranking','Characters','Revisions','Hours
'),1:70))
> 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
Ranking Characters Revisions Hours\r
1 1 21221 6086 7
2 2 7106 2273 6
3 3 4794 2601 6
4 4 6776 2708 6
5 5 7122 3321 6
6 6 12025 2106 6
7 7 18945 6353 5
8 8 24380 7912 5
9 9 6942 1905 5
10 10 7370 3841 5
11 11 13143 2737 5
12 12 11831 4449 5
13 13 3102 1253 5
14 14 7355 2372 5
15 15 7749 3303 5
16 16 10306 2876 5
17 17 7860 2939 4
18 18 9138 2149 4
19 19 6545 2545 4
20 20 9829 3734 4
21 21 9460 1911 4
22 22 30631 875 4
23 23 4807 1466 4
24 24 4618 1562 4
25 25 7304 1203 4
26 26 12621 2860 4
27 27 9320 3245 4
28 28 9569 2668 4
29 29 8480 2111 4
30 30 4718 942 4
31 31 7342 2690 4
32 32 8230 3675 4
33 33 10007 3798 4
34 34 9176 672 4
35 35 2334 523 4
36 36 5459 2352 3
37 37 5440 1476 3
38 38 9027 2646 3
39 39 9852 2500 3
40 40 3122 720 3
41 41 8863 2483 3
42 42 4286 1550 3
43 43 9119 2551 3
44 44 7297 2544 3
45 45 8009 2870 3
46 46 5049 1288 3
47 47 5057 1474 3
48 48 4372 1521 3
49 49 3393 1331 3
50 50 7586 1685 3
51 51 7440 2755 3
52 52 8636 1268 3
53 53 10846 1266 3
54 54 6843 1564 3
55 55 5102 1062 3
56 56 4183 1306 3
57 57 10083 4566 3
58 58 5273 1655 3
59 59 6433 1423 3
60 60 6116 474 3
61 61 5887 725 3
62 62 7415 1061 2
63 63 7953 687 2
64 64 8447 2411 2
65 65 4733 1705 2
66 66 11478 2422 2
67 67 4234 1900 2
68 68 11089 2848 2
69 69 6390 1398 2
70 70 7737 2625 2
> k <- length(x[1,])
> df <- as.data.frame(x)
> (mylm <- lm(df))
Call:
lm(formula = df)
Coefficients:
(Intercept) Characters Revisions `Hours\\r`
9.649e+01 -1.141e-04 -3.876e-04 -1.598e+01
> (mysum <- summary(mylm))
Call:
lm(formula = df)
Residuals:
Min 1Q Median 3Q Max
-13.5173 -4.7695 -0.3675 4.2772 21.1773
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 9.649e+01 2.921e+00 33.037 <2e-16 ***
Characters -1.141e-04 2.242e-04 -0.509 0.613
Revisions -3.876e-04 8.230e-04 -0.471 0.639
`Hours\\r` -1.598e+01 8.145e-01 -19.623 <2e-16 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 7.077 on 66 degrees of freedom
Multiple R-squared: 0.8843, Adjusted R-squared: 0.8791
F-statistic: 168.2 on 3 and 66 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,] 2.275577e-02 4.551154e-02 0.977244231
[2,] 4.970907e-03 9.941814e-03 0.995029093
[3,] 1.244508e-03 2.489016e-03 0.998755492
[4,] 1.034344e-03 2.068688e-03 0.998965656
[5,] 3.527939e-04 7.055878e-04 0.999647206
[6,] 4.800418e-04 9.600836e-04 0.999519958
[7,] 4.074384e-04 8.148768e-04 0.999592562
[8,] 4.400187e-04 8.800375e-04 0.999559981
[9,] 8.518153e-04 1.703631e-03 0.999148185
[10,] 1.370307e-03 2.740613e-03 0.998629693
[11,] 6.898755e-04 1.379751e-03 0.999310124
[12,] 3.446193e-04 6.892386e-04 0.999655381
[13,] 1.883773e-04 3.767545e-04 0.999811623
[14,] 1.265984e-04 2.531968e-04 0.999873402
[15,] 7.300153e-05 1.460031e-04 0.999926998
[16,] 3.284060e-05 6.568119e-05 0.999967159
[17,] 3.730359e-05 7.460718e-05 0.999962696
[18,] 5.154466e-05 1.030893e-04 0.999948455
[19,] 6.933661e-05 1.386732e-04 0.999930663
[20,] 1.665333e-04 3.330667e-04 0.999833467
[21,] 4.697770e-04 9.395540e-04 0.999530223
[22,] 1.006258e-03 2.012516e-03 0.998993742
[23,] 1.889053e-03 3.778106e-03 0.998110947
[24,] 2.850131e-03 5.700262e-03 0.997149869
[25,] 5.856680e-03 1.171336e-02 0.994143320
[26,] 1.257658e-02 2.515316e-02 0.987423418
[27,] 2.394821e-02 4.789641e-02 0.976051795
[28,] 3.355121e-02 6.710242e-02 0.966448788
[29,] 4.218799e-02 8.437598e-02 0.957812012
[30,] 4.941092e-02 9.882184e-02 0.950589081
[31,] 6.904510e-02 1.380902e-01 0.930954896
[32,] 8.868656e-02 1.773731e-01 0.911313438
[33,] 1.273517e-01 2.547033e-01 0.872648340
[34,] 1.828943e-01 3.657886e-01 0.817105677
[35,] 2.544048e-01 5.088096e-01 0.745595202
[36,] 3.426536e-01 6.853071e-01 0.657346428
[37,] 4.564098e-01 9.128196e-01 0.543590193
[38,] 5.753945e-01 8.492110e-01 0.424605504
[39,] 7.099182e-01 5.801635e-01 0.290081761
[40,] 8.034370e-01 3.931261e-01 0.196563046
[41,] 8.787510e-01 2.424979e-01 0.121248967
[42,] 9.333189e-01 1.333623e-01 0.066681126
[43,] 9.728300e-01 5.433999e-02 0.027169996
[44,] 9.864386e-01 2.712281e-02 0.013561407
[45,] 9.951480e-01 9.704045e-03 0.004852023
[46,] 9.965314e-01 6.937279e-03 0.003468639
[47,] 9.968986e-01 6.202886e-03 0.003101443
[48,] 9.976235e-01 4.753007e-03 0.002376504
[49,] 9.976616e-01 4.676809e-03 0.002338404
[50,] 9.974023e-01 5.195433e-03 0.002597717
[51,] 9.983774e-01 3.245108e-03 0.001622554
[52,] 9.983445e-01 3.311068e-03 0.001655534
[53,] 9.978446e-01 4.310858e-03 0.002155429
[54,] 9.942542e-01 1.149153e-02 0.005745766
[55,] 9.847022e-01 3.059562e-02 0.015297810
[56,] 9.655083e-01 6.898344e-02 0.034491722
[57,] 9.112089e-01 1.775821e-01 0.088791060
> postscript(file="/var/fisher/rcomp/tmp/1j19m1352156420.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/fisher/rcomp/tmp/2s0fi1352156420.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/fisher/rcomp/tmp/31ser1352156420.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/fisher/rcomp/tmp/432xn1352156420.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/fisher/rcomp/tmp/50gjg1352156420.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 = 70
Frequency = 1
1 2 3 4 5 6
21.1772877 3.1055307 3.9688639 5.2364707 6.5135259 7.6020424
7 8 9 10 11 12
-4.9459317 -2.7216070 -6.0393147 -4.2401531 -3.0093537 -1.4955330
13 14 15 16 17 18
-2.7301344 -0.8111997 0.5945783 1.7208294 -13.5173406 -12.6777042
19 20 21 22 23 24
-11.8200772 -9.9845714 -9.7332063 -6.7192113 -8.4365593 -7.4209170
25 26 27 28 29 30
-6.2535931 -4.0047501 -3.2321659 -2.4273821 -1.7675065 -1.6497987
31 32 33 34 35 36
0.3270539 1.8101237 3.0605419 2.7541957 2.9158075 -11.0022945
37 38 39 40 41 42
-10.3439705 -8.4812575 -7.4437135 -7.9014440 -5.5631425 -5.4469567
43 44 45 46 47 48
-3.5075795 -2.7181743 -1.5105917 -1.4614443 -0.3884442 0.5516161
49 50 51 52 53 54
1.3662790 2.9818796 4.3799178 4.9400646 6.1914403 6.8502112
55 56 57 58 59 60
7.4570127 8.4467252 11.3833550 10.7063499 11.7487852 12.3448163
61 62 63 64 65 66
13.4159677 -1.2629815 -0.3465480 1.3779798 1.6806075 3.7280662
67 68 69 70
3.6992494 5.8487864 5.7506807 7.3799113
> postscript(file="/var/fisher/rcomp/tmp/6w3nc1352156420.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 = 70
Frequency = 1
lag(myerror, k = 1) myerror
0 21.1772877 NA
1 3.1055307 21.1772877
2 3.9688639 3.1055307
3 5.2364707 3.9688639
4 6.5135259 5.2364707
5 7.6020424 6.5135259
6 -4.9459317 7.6020424
7 -2.7216070 -4.9459317
8 -6.0393147 -2.7216070
9 -4.2401531 -6.0393147
10 -3.0093537 -4.2401531
11 -1.4955330 -3.0093537
12 -2.7301344 -1.4955330
13 -0.8111997 -2.7301344
14 0.5945783 -0.8111997
15 1.7208294 0.5945783
16 -13.5173406 1.7208294
17 -12.6777042 -13.5173406
18 -11.8200772 -12.6777042
19 -9.9845714 -11.8200772
20 -9.7332063 -9.9845714
21 -6.7192113 -9.7332063
22 -8.4365593 -6.7192113
23 -7.4209170 -8.4365593
24 -6.2535931 -7.4209170
25 -4.0047501 -6.2535931
26 -3.2321659 -4.0047501
27 -2.4273821 -3.2321659
28 -1.7675065 -2.4273821
29 -1.6497987 -1.7675065
30 0.3270539 -1.6497987
31 1.8101237 0.3270539
32 3.0605419 1.8101237
33 2.7541957 3.0605419
34 2.9158075 2.7541957
35 -11.0022945 2.9158075
36 -10.3439705 -11.0022945
37 -8.4812575 -10.3439705
38 -7.4437135 -8.4812575
39 -7.9014440 -7.4437135
40 -5.5631425 -7.9014440
41 -5.4469567 -5.5631425
42 -3.5075795 -5.4469567
43 -2.7181743 -3.5075795
44 -1.5105917 -2.7181743
45 -1.4614443 -1.5105917
46 -0.3884442 -1.4614443
47 0.5516161 -0.3884442
48 1.3662790 0.5516161
49 2.9818796 1.3662790
50 4.3799178 2.9818796
51 4.9400646 4.3799178
52 6.1914403 4.9400646
53 6.8502112 6.1914403
54 7.4570127 6.8502112
55 8.4467252 7.4570127
56 11.3833550 8.4467252
57 10.7063499 11.3833550
58 11.7487852 10.7063499
59 12.3448163 11.7487852
60 13.4159677 12.3448163
61 -1.2629815 13.4159677
62 -0.3465480 -1.2629815
63 1.3779798 -0.3465480
64 1.6806075 1.3779798
65 3.7280662 1.6806075
66 3.6992494 3.7280662
67 5.8487864 3.6992494
68 5.7506807 5.8487864
69 7.3799113 5.7506807
70 NA 7.3799113
> dum1 <- dum[2:length(myerror),]
> dum1
lag(myerror, k = 1) myerror
[1,] 3.1055307 21.1772877
[2,] 3.9688639 3.1055307
[3,] 5.2364707 3.9688639
[4,] 6.5135259 5.2364707
[5,] 7.6020424 6.5135259
[6,] -4.9459317 7.6020424
[7,] -2.7216070 -4.9459317
[8,] -6.0393147 -2.7216070
[9,] -4.2401531 -6.0393147
[10,] -3.0093537 -4.2401531
[11,] -1.4955330 -3.0093537
[12,] -2.7301344 -1.4955330
[13,] -0.8111997 -2.7301344
[14,] 0.5945783 -0.8111997
[15,] 1.7208294 0.5945783
[16,] -13.5173406 1.7208294
[17,] -12.6777042 -13.5173406
[18,] -11.8200772 -12.6777042
[19,] -9.9845714 -11.8200772
[20,] -9.7332063 -9.9845714
[21,] -6.7192113 -9.7332063
[22,] -8.4365593 -6.7192113
[23,] -7.4209170 -8.4365593
[24,] -6.2535931 -7.4209170
[25,] -4.0047501 -6.2535931
[26,] -3.2321659 -4.0047501
[27,] -2.4273821 -3.2321659
[28,] -1.7675065 -2.4273821
[29,] -1.6497987 -1.7675065
[30,] 0.3270539 -1.6497987
[31,] 1.8101237 0.3270539
[32,] 3.0605419 1.8101237
[33,] 2.7541957 3.0605419
[34,] 2.9158075 2.7541957
[35,] -11.0022945 2.9158075
[36,] -10.3439705 -11.0022945
[37,] -8.4812575 -10.3439705
[38,] -7.4437135 -8.4812575
[39,] -7.9014440 -7.4437135
[40,] -5.5631425 -7.9014440
[41,] -5.4469567 -5.5631425
[42,] -3.5075795 -5.4469567
[43,] -2.7181743 -3.5075795
[44,] -1.5105917 -2.7181743
[45,] -1.4614443 -1.5105917
[46,] -0.3884442 -1.4614443
[47,] 0.5516161 -0.3884442
[48,] 1.3662790 0.5516161
[49,] 2.9818796 1.3662790
[50,] 4.3799178 2.9818796
[51,] 4.9400646 4.3799178
[52,] 6.1914403 4.9400646
[53,] 6.8502112 6.1914403
[54,] 7.4570127 6.8502112
[55,] 8.4467252 7.4570127
[56,] 11.3833550 8.4467252
[57,] 10.7063499 11.3833550
[58,] 11.7487852 10.7063499
[59,] 12.3448163 11.7487852
[60,] 13.4159677 12.3448163
[61,] -1.2629815 13.4159677
[62,] -0.3465480 -1.2629815
[63,] 1.3779798 -0.3465480
[64,] 1.6806075 1.3779798
[65,] 3.7280662 1.6806075
[66,] 3.6992494 3.7280662
[67,] 5.8487864 3.6992494
[68,] 5.7506807 5.8487864
[69,] 7.3799113 5.7506807
> z <- as.data.frame(dum1)
> z
lag(myerror, k = 1) myerror
1 3.1055307 21.1772877
2 3.9688639 3.1055307
3 5.2364707 3.9688639
4 6.5135259 5.2364707
5 7.6020424 6.5135259
6 -4.9459317 7.6020424
7 -2.7216070 -4.9459317
8 -6.0393147 -2.7216070
9 -4.2401531 -6.0393147
10 -3.0093537 -4.2401531
11 -1.4955330 -3.0093537
12 -2.7301344 -1.4955330
13 -0.8111997 -2.7301344
14 0.5945783 -0.8111997
15 1.7208294 0.5945783
16 -13.5173406 1.7208294
17 -12.6777042 -13.5173406
18 -11.8200772 -12.6777042
19 -9.9845714 -11.8200772
20 -9.7332063 -9.9845714
21 -6.7192113 -9.7332063
22 -8.4365593 -6.7192113
23 -7.4209170 -8.4365593
24 -6.2535931 -7.4209170
25 -4.0047501 -6.2535931
26 -3.2321659 -4.0047501
27 -2.4273821 -3.2321659
28 -1.7675065 -2.4273821
29 -1.6497987 -1.7675065
30 0.3270539 -1.6497987
31 1.8101237 0.3270539
32 3.0605419 1.8101237
33 2.7541957 3.0605419
34 2.9158075 2.7541957
35 -11.0022945 2.9158075
36 -10.3439705 -11.0022945
37 -8.4812575 -10.3439705
38 -7.4437135 -8.4812575
39 -7.9014440 -7.4437135
40 -5.5631425 -7.9014440
41 -5.4469567 -5.5631425
42 -3.5075795 -5.4469567
43 -2.7181743 -3.5075795
44 -1.5105917 -2.7181743
45 -1.4614443 -1.5105917
46 -0.3884442 -1.4614443
47 0.5516161 -0.3884442
48 1.3662790 0.5516161
49 2.9818796 1.3662790
50 4.3799178 2.9818796
51 4.9400646 4.3799178
52 6.1914403 4.9400646
53 6.8502112 6.1914403
54 7.4570127 6.8502112
55 8.4467252 7.4570127
56 11.3833550 8.4467252
57 10.7063499 11.3833550
58 11.7487852 10.7063499
59 12.3448163 11.7487852
60 13.4159677 12.3448163
61 -1.2629815 13.4159677
62 -0.3465480 -1.2629815
63 1.3779798 -0.3465480
64 1.6806075 1.3779798
65 3.7280662 1.6806075
66 3.6992494 3.7280662
67 5.8487864 3.6992494
68 5.7506807 5.8487864
69 7.3799113 5.7506807
> 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/fisher/rcomp/tmp/70w8g1352156420.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/fisher/rcomp/tmp/8o2ls1352156420.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/fisher/rcomp/tmp/9l8fr1352156420.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/fisher/rcomp/tmp/104r3f1352156420.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/fisher/rcomp/createtable file can be downloaded at http://www.wessa.net/cretab
> load(file="/var/fisher/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/fisher/rcomp/tmp/11wlww1352156420.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/fisher/rcomp/tmp/12x2l31352156420.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/fisher/rcomp/tmp/133qlv1352156421.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/fisher/rcomp/tmp/14dcdx1352156421.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/fisher/rcomp/tmp/150zpm1352156421.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/fisher/rcomp/tmp/16l1hp1352156421.tab")
+ }
>
> try(system("convert tmp/1j19m1352156420.ps tmp/1j19m1352156420.png",intern=TRUE))
character(0)
> try(system("convert tmp/2s0fi1352156420.ps tmp/2s0fi1352156420.png",intern=TRUE))
character(0)
> try(system("convert tmp/31ser1352156420.ps tmp/31ser1352156420.png",intern=TRUE))
character(0)
> try(system("convert tmp/432xn1352156420.ps tmp/432xn1352156420.png",intern=TRUE))
character(0)
> try(system("convert tmp/50gjg1352156420.ps tmp/50gjg1352156420.png",intern=TRUE))
character(0)
> try(system("convert tmp/6w3nc1352156420.ps tmp/6w3nc1352156420.png",intern=TRUE))
character(0)
> try(system("convert tmp/70w8g1352156420.ps tmp/70w8g1352156420.png",intern=TRUE))
character(0)
> try(system("convert tmp/8o2ls1352156420.ps tmp/8o2ls1352156420.png",intern=TRUE))
character(0)
> try(system("convert tmp/9l8fr1352156420.ps tmp/9l8fr1352156420.png",intern=TRUE))
character(0)
> try(system("convert tmp/104r3f1352156420.ps tmp/104r3f1352156420.png",intern=TRUE))
character(0)
>
>
> proc.time()
user system elapsed
6.405 1.201 7.606