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(2649.2
+ ,31077
+ ,0
+ ,2579.4
+ ,31293
+ ,0
+ ,2504.6
+ ,30236
+ ,0
+ ,2462.3
+ ,30160
+ ,0
+ ,2467.4
+ ,32436
+ ,0
+ ,2446.7
+ ,30695
+ ,0
+ ,2656.3
+ ,27525
+ ,0
+ ,2626.2
+ ,26434
+ ,0
+ ,2482.6
+ ,25739
+ ,0
+ ,2539.9
+ ,25204
+ ,0
+ ,2502.7
+ ,24977
+ ,0
+ ,2466.9
+ ,24320
+ ,0
+ ,2513.2
+ ,22680
+ ,1
+ ,2443.3
+ ,22052
+ ,1
+ ,2293.4
+ ,21467
+ ,1
+ ,2070.8
+ ,21383
+ ,1
+ ,2029.6
+ ,21777
+ ,1
+ ,2052
+ ,21928
+ ,1
+ ,1864.4
+ ,21814
+ ,1
+ ,1670.1
+ ,22937
+ ,1
+ ,1811
+ ,23595
+ ,1
+ ,1905.4
+ ,20830
+ ,1
+ ,1862.8
+ ,19650
+ ,1
+ ,2014.5
+ ,19195
+ ,1
+ ,2197.8
+ ,19644
+ ,1
+ ,2962.3
+ ,18483
+ ,0
+ ,3047
+ ,18079
+ ,0
+ ,3032.6
+ ,19178
+ ,0
+ ,3504.4
+ ,18391
+ ,0
+ ,3801.1
+ ,18441
+ ,0
+ ,3857.6
+ ,18584
+ ,0
+ ,3674.4
+ ,20108
+ ,0
+ ,3721
+ ,20148
+ ,0
+ ,3844.5
+ ,19394
+ ,0
+ ,4116.7
+ ,17745
+ ,0
+ ,4105.2
+ ,17696
+ ,0
+ ,4435.2
+ ,17032
+ ,0
+ ,4296.5
+ ,16438
+ ,0
+ ,4202.5
+ ,15683
+ ,0
+ ,4562.8
+ ,15594
+ ,0
+ ,4621.4
+ ,15713
+ ,0
+ ,4697
+ ,15937
+ ,0
+ ,4591.3
+ ,16171
+ ,0
+ ,4357
+ ,15928
+ ,0
+ ,4502.6
+ ,16348
+ ,0
+ ,4443.9
+ ,15579
+ ,0
+ ,4290.9
+ ,15305
+ ,0
+ ,4199.8
+ ,15648
+ ,0
+ ,4138.5
+ ,14954
+ ,0
+ ,3970.1
+ ,15137
+ ,0
+ ,3862.3
+ ,15839
+ ,0
+ ,3701.6
+ ,16050
+ ,0
+ ,3570.12
+ ,15168
+ ,0
+ ,3801.06
+ ,17064
+ ,0
+ ,3895.51
+ ,16005
+ ,0
+ ,3917.96
+ ,14886
+ ,0
+ ,3813.06
+ ,14931
+ ,0
+ ,3667.03
+ ,14544
+ ,0
+ ,3494.17
+ ,13812
+ ,0
+ ,3364
+ ,13031
+ ,0
+ ,3295.3
+ ,12574
+ ,0
+ ,3277.0
+ ,11964
+ ,0
+ ,3257.2
+ ,11451
+ ,0
+ ,3161.7
+ ,11346
+ ,0
+ ,3097.3
+ ,11353
+ ,0
+ ,3061.3
+ ,10702
+ ,0
+ ,3119.3
+ ,10646
+ ,0
+ ,3106.22
+ ,10556
+ ,0
+ ,3080.58
+ ,10463
+ ,0
+ ,2981.85
+ ,10407
+ ,0)
+ ,dim=c(3
+ ,70)
+ ,dimnames=list(c('Bel20'
+ ,'Goudprijs'
+ ,'Crisis')
+ ,1:70))
> y <- array(NA,dim=c(3,70),dimnames=list(c('Bel20','Goudprijs','Crisis'),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 = 'Include Monthly 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.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
Bel20 Goudprijs Crisis M1 M2 M3 M4 M5 M6 M7 M8 M9 M10 M11
1 2649.20 31077 0 1 0 0 0 0 0 0 0 0 0 0
2 2579.40 31293 0 0 1 0 0 0 0 0 0 0 0 0
3 2504.60 30236 0 0 0 1 0 0 0 0 0 0 0 0
4 2462.30 30160 0 0 0 0 1 0 0 0 0 0 0 0
5 2467.40 32436 0 0 0 0 0 1 0 0 0 0 0 0
6 2446.70 30695 0 0 0 0 0 0 1 0 0 0 0 0
7 2656.30 27525 0 0 0 0 0 0 0 1 0 0 0 0
8 2626.20 26434 0 0 0 0 0 0 0 0 1 0 0 0
9 2482.60 25739 0 0 0 0 0 0 0 0 0 1 0 0
10 2539.90 25204 0 0 0 0 0 0 0 0 0 0 1 0
11 2502.70 24977 0 0 0 0 0 0 0 0 0 0 0 1
12 2466.90 24320 0 0 0 0 0 0 0 0 0 0 0 0
13 2513.20 22680 1 1 0 0 0 0 0 0 0 0 0 0
14 2443.30 22052 1 0 1 0 0 0 0 0 0 0 0 0
15 2293.40 21467 1 0 0 1 0 0 0 0 0 0 0 0
16 2070.80 21383 1 0 0 0 1 0 0 0 0 0 0 0
17 2029.60 21777 1 0 0 0 0 1 0 0 0 0 0 0
18 2052.00 21928 1 0 0 0 0 0 1 0 0 0 0 0
19 1864.40 21814 1 0 0 0 0 0 0 1 0 0 0 0
20 1670.10 22937 1 0 0 0 0 0 0 0 1 0 0 0
21 1811.00 23595 1 0 0 0 0 0 0 0 0 1 0 0
22 1905.40 20830 1 0 0 0 0 0 0 0 0 0 1 0
23 1862.80 19650 1 0 0 0 0 0 0 0 0 0 0 1
24 2014.50 19195 1 0 0 0 0 0 0 0 0 0 0 0
25 2197.80 19644 1 1 0 0 0 0 0 0 0 0 0 0
26 2962.30 18483 0 0 1 0 0 0 0 0 0 0 0 0
27 3047.00 18079 0 0 0 1 0 0 0 0 0 0 0 0
28 3032.60 19178 0 0 0 0 1 0 0 0 0 0 0 0
29 3504.40 18391 0 0 0 0 0 1 0 0 0 0 0 0
30 3801.10 18441 0 0 0 0 0 0 1 0 0 0 0 0
31 3857.60 18584 0 0 0 0 0 0 0 1 0 0 0 0
32 3674.40 20108 0 0 0 0 0 0 0 0 1 0 0 0
33 3721.00 20148 0 0 0 0 0 0 0 0 0 1 0 0
34 3844.50 19394 0 0 0 0 0 0 0 0 0 0 1 0
35 4116.70 17745 0 0 0 0 0 0 0 0 0 0 0 1
36 4105.20 17696 0 0 0 0 0 0 0 0 0 0 0 0
37 4435.20 17032 0 1 0 0 0 0 0 0 0 0 0 0
38 4296.50 16438 0 0 1 0 0 0 0 0 0 0 0 0
39 4202.50 15683 0 0 0 1 0 0 0 0 0 0 0 0
40 4562.80 15594 0 0 0 0 1 0 0 0 0 0 0 0
41 4621.40 15713 0 0 0 0 0 1 0 0 0 0 0 0
42 4697.00 15937 0 0 0 0 0 0 1 0 0 0 0 0
43 4591.30 16171 0 0 0 0 0 0 0 1 0 0 0 0
44 4357.00 15928 0 0 0 0 0 0 0 0 1 0 0 0
45 4502.60 16348 0 0 0 0 0 0 0 0 0 1 0 0
46 4443.90 15579 0 0 0 0 0 0 0 0 0 0 1 0
47 4290.90 15305 0 0 0 0 0 0 0 0 0 0 0 1
48 4199.80 15648 0 0 0 0 0 0 0 0 0 0 0 0
49 4138.50 14954 0 1 0 0 0 0 0 0 0 0 0 0
50 3970.10 15137 0 0 1 0 0 0 0 0 0 0 0 0
51 3862.30 15839 0 0 0 1 0 0 0 0 0 0 0 0
52 3701.60 16050 0 0 0 0 1 0 0 0 0 0 0 0
53 3570.12 15168 0 0 0 0 0 1 0 0 0 0 0 0
54 3801.06 17064 0 0 0 0 0 0 1 0 0 0 0 0
55 3895.51 16005 0 0 0 0 0 0 0 1 0 0 0 0
56 3917.96 14886 0 0 0 0 0 0 0 0 1 0 0 0
57 3813.06 14931 0 0 0 0 0 0 0 0 0 1 0 0
58 3667.03 14544 0 0 0 0 0 0 0 0 0 0 1 0
59 3494.17 13812 0 0 0 0 0 0 0 0 0 0 0 1
60 3364.00 13031 0 0 0 0 0 0 0 0 0 0 0 0
61 3295.30 12574 0 1 0 0 0 0 0 0 0 0 0 0
62 3277.00 11964 0 0 1 0 0 0 0 0 0 0 0 0
63 3257.20 11451 0 0 0 1 0 0 0 0 0 0 0 0
64 3161.70 11346 0 0 0 0 1 0 0 0 0 0 0 0
65 3097.30 11353 0 0 0 0 0 1 0 0 0 0 0 0
66 3061.30 10702 0 0 0 0 0 0 1 0 0 0 0 0
67 3119.30 10646 0 0 0 0 0 0 0 1 0 0 0 0
68 3106.22 10556 0 0 0 0 0 0 0 0 1 0 0 0
69 3080.58 10463 0 0 0 0 0 0 0 0 0 1 0 0
70 2981.85 10407 0 0 0 0 0 0 0 0 0 0 1 0
> k <- length(x[1,])
> df <- as.data.frame(x)
> (mylm <- lm(df))
Call:
lm(formula = df)
Coefficients:
(Intercept) Goudprijs Crisis M1 M2 M3
4.594e+03 -6.184e-02 -1.262e+03 2.470e+02 5.992e+01 -2.727e+01
M4 M5 M6 M7 M8 M9
-4.661e+01 1.474e+01 1.088e+02 8.825e+01 -1.610e+01 -2.405e+00
M10 M11
-6.139e+01 4.315e+01
> (mysum <- summary(mylm))
Call:
lm(formula = df)
Residuals:
Min 1Q Median 3Q Max
-979.85 -350.85 -46.72 338.89 984.20
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 4.594e+03 3.476e+02 13.217 < 2e-16 ***
Goudprijs -6.184e-02 1.301e-02 -4.754 1.44e-05 ***
Crisis -1.262e+03 1.868e+02 -6.755 8.80e-09 ***
M1 2.470e+02 3.564e+02 0.693 0.491
M2 5.992e+01 3.557e+02 0.168 0.867
M3 -2.727e+01 3.555e+02 -0.077 0.939
M4 -4.661e+01 3.555e+02 -0.131 0.896
M5 1.474e+01 3.556e+02 0.041 0.967
M6 1.088e+02 3.556e+02 0.306 0.761
M7 8.825e+01 3.553e+02 0.248 0.805
M8 -1.610e+01 3.553e+02 -0.045 0.964
M9 -2.405e+00 3.554e+02 -0.007 0.995
M10 -6.139e+01 3.553e+02 -0.173 0.863
M11 4.315e+01 3.710e+02 0.116 0.908
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Residual standard error: 586.6 on 56 degrees of freedom
Multiple R-squared: 0.6123, Adjusted R-squared: 0.5223
F-statistic: 6.802 on 13 and 56 DF, p-value: 1.320e-07
> 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,] 8.326970e-03 1.665394e-02 0.9916730296
[2,] 2.629341e-03 5.258682e-03 0.9973706591
[3,] 5.075126e-03 1.015025e-02 0.9949248740
[4,] 1.349320e-03 2.698640e-03 0.9986506802
[5,] 8.232962e-04 1.646592e-03 0.9991767038
[6,] 2.057133e-04 4.114265e-04 0.9997942867
[7,] 4.827436e-05 9.654871e-05 0.9999517256
[8,] 1.554944e-05 3.109889e-05 0.9999844506
[9,] 3.467255e-05 6.934510e-05 0.9999653275
[10,] 3.659998e-05 7.319997e-05 0.9999634000
[11,] 1.994277e-05 3.988553e-05 0.9999800572
[12,] 2.680358e-05 5.360716e-05 0.9999731964
[13,] 1.818609e-04 3.637218e-04 0.9998181391
[14,] 1.975185e-03 3.950370e-03 0.9980248152
[15,] 9.736752e-03 1.947350e-02 0.9902632478
[16,] 4.295706e-02 8.591413e-02 0.9570429352
[17,] 1.670392e-01 3.340784e-01 0.8329607882
[18,] 4.348909e-01 8.697818e-01 0.5651090918
[19,] 6.120389e-01 7.759221e-01 0.3879610642
[20,] 7.064732e-01 5.870537e-01 0.2935268472
[21,] 6.544444e-01 6.911112e-01 0.3455555769
[22,] 6.024958e-01 7.950084e-01 0.3975042103
[23,] 5.440719e-01 9.118562e-01 0.4559281000
[24,] 6.761026e-01 6.477948e-01 0.3238973791
[25,] 7.930440e-01 4.139121e-01 0.2069560467
[26,] 9.255236e-01 1.489528e-01 0.0744764065
[27,] 9.523384e-01 9.532325e-02 0.0476616270
[28,] 9.370385e-01 1.259231e-01 0.0629615264
[29,] 9.430708e-01 1.138583e-01 0.0569291588
[30,] 9.779245e-01 4.415093e-02 0.0220754674
[31,] 9.898728e-01 2.025437e-02 0.0101271838
[32,] 9.934382e-01 1.312367e-02 0.0065618364
[33,] 9.991133e-01 1.773498e-03 0.0008867489
[34,] 9.990554e-01 1.889175e-03 0.0009445874
[35,] 9.959804e-01 8.039184e-03 0.0040195922
[36,] 9.891892e-01 2.162165e-02 0.0108108255
[37,] 9.668537e-01 6.629268e-02 0.0331463396
> postscript(file="/var/www/html/freestat/rcomp/tmp/185jx1291107212.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/285jx1291107212.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/3jei01291107212.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/4jei01291107212.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/5jei01291107212.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 7
-270.21392 -139.55831 -192.53361 -220.18590 -135.69672 -358.14644 -323.99318
8 9 10 11 12 13 14
-317.20746 -517.47561 -434.27615 -590.05029 -623.32771 336.12926 414.49450
15 16 17 18 19 20 21
315.60624 107.15925 28.97130 -33.38292 -207.45597 -227.96325 -60.06602
22 23 24 25 26 27 28
-77.66295 -297.76769 -131.05404 -167.00770 -548.78955 -401.88533 -328.97921
29 30 31 32 33 34 35
32.80341 238.50366 324.42271 339.81203 375.19404 511.05121 576.74494
36 37 38 39 40 41 42
605.36438 647.28621 658.95390 605.45336 979.59718 984.20407 979.56395
43 44 45 46 47 48 49
908.91016 763.93361 921.81365 874.54327 600.06280 573.32232 222.08904
50 51 52 53 54 55 56
252.10404 274.89992 146.59483 -100.77706 153.31419 202.85523 260.45952
57 58 59 60 61 62 63
144.65070 33.67203 -288.98976 -424.30496 -768.28289 -637.20458 -601.54059
64 65 66 67 68 69 70
-684.18615 -809.50500 -979.85244 -904.73895 -819.03445 -864.11676 -907.32740
> postscript(file="/var/www/html/freestat/rcomp/tmp/6u6il1291107212.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 -270.21392 NA
1 -139.55831 -270.21392
2 -192.53361 -139.55831
3 -220.18590 -192.53361
4 -135.69672 -220.18590
5 -358.14644 -135.69672
6 -323.99318 -358.14644
7 -317.20746 -323.99318
8 -517.47561 -317.20746
9 -434.27615 -517.47561
10 -590.05029 -434.27615
11 -623.32771 -590.05029
12 336.12926 -623.32771
13 414.49450 336.12926
14 315.60624 414.49450
15 107.15925 315.60624
16 28.97130 107.15925
17 -33.38292 28.97130
18 -207.45597 -33.38292
19 -227.96325 -207.45597
20 -60.06602 -227.96325
21 -77.66295 -60.06602
22 -297.76769 -77.66295
23 -131.05404 -297.76769
24 -167.00770 -131.05404
25 -548.78955 -167.00770
26 -401.88533 -548.78955
27 -328.97921 -401.88533
28 32.80341 -328.97921
29 238.50366 32.80341
30 324.42271 238.50366
31 339.81203 324.42271
32 375.19404 339.81203
33 511.05121 375.19404
34 576.74494 511.05121
35 605.36438 576.74494
36 647.28621 605.36438
37 658.95390 647.28621
38 605.45336 658.95390
39 979.59718 605.45336
40 984.20407 979.59718
41 979.56395 984.20407
42 908.91016 979.56395
43 763.93361 908.91016
44 921.81365 763.93361
45 874.54327 921.81365
46 600.06280 874.54327
47 573.32232 600.06280
48 222.08904 573.32232
49 252.10404 222.08904
50 274.89992 252.10404
51 146.59483 274.89992
52 -100.77706 146.59483
53 153.31419 -100.77706
54 202.85523 153.31419
55 260.45952 202.85523
56 144.65070 260.45952
57 33.67203 144.65070
58 -288.98976 33.67203
59 -424.30496 -288.98976
60 -768.28289 -424.30496
61 -637.20458 -768.28289
62 -601.54059 -637.20458
63 -684.18615 -601.54059
64 -809.50500 -684.18615
65 -979.85244 -809.50500
66 -904.73895 -979.85244
67 -819.03445 -904.73895
68 -864.11676 -819.03445
69 -907.32740 -864.11676
70 NA -907.32740
> dum1 <- dum[2:length(myerror),]
> dum1
lag(myerror, k = 1) myerror
[1,] -139.55831 -270.21392
[2,] -192.53361 -139.55831
[3,] -220.18590 -192.53361
[4,] -135.69672 -220.18590
[5,] -358.14644 -135.69672
[6,] -323.99318 -358.14644
[7,] -317.20746 -323.99318
[8,] -517.47561 -317.20746
[9,] -434.27615 -517.47561
[10,] -590.05029 -434.27615
[11,] -623.32771 -590.05029
[12,] 336.12926 -623.32771
[13,] 414.49450 336.12926
[14,] 315.60624 414.49450
[15,] 107.15925 315.60624
[16,] 28.97130 107.15925
[17,] -33.38292 28.97130
[18,] -207.45597 -33.38292
[19,] -227.96325 -207.45597
[20,] -60.06602 -227.96325
[21,] -77.66295 -60.06602
[22,] -297.76769 -77.66295
[23,] -131.05404 -297.76769
[24,] -167.00770 -131.05404
[25,] -548.78955 -167.00770
[26,] -401.88533 -548.78955
[27,] -328.97921 -401.88533
[28,] 32.80341 -328.97921
[29,] 238.50366 32.80341
[30,] 324.42271 238.50366
[31,] 339.81203 324.42271
[32,] 375.19404 339.81203
[33,] 511.05121 375.19404
[34,] 576.74494 511.05121
[35,] 605.36438 576.74494
[36,] 647.28621 605.36438
[37,] 658.95390 647.28621
[38,] 605.45336 658.95390
[39,] 979.59718 605.45336
[40,] 984.20407 979.59718
[41,] 979.56395 984.20407
[42,] 908.91016 979.56395
[43,] 763.93361 908.91016
[44,] 921.81365 763.93361
[45,] 874.54327 921.81365
[46,] 600.06280 874.54327
[47,] 573.32232 600.06280
[48,] 222.08904 573.32232
[49,] 252.10404 222.08904
[50,] 274.89992 252.10404
[51,] 146.59483 274.89992
[52,] -100.77706 146.59483
[53,] 153.31419 -100.77706
[54,] 202.85523 153.31419
[55,] 260.45952 202.85523
[56,] 144.65070 260.45952
[57,] 33.67203 144.65070
[58,] -288.98976 33.67203
[59,] -424.30496 -288.98976
[60,] -768.28289 -424.30496
[61,] -637.20458 -768.28289
[62,] -601.54059 -637.20458
[63,] -684.18615 -601.54059
[64,] -809.50500 -684.18615
[65,] -979.85244 -809.50500
[66,] -904.73895 -979.85244
[67,] -819.03445 -904.73895
[68,] -864.11676 -819.03445
[69,] -907.32740 -864.11676
> z <- as.data.frame(dum1)
> z
lag(myerror, k = 1) myerror
1 -139.55831 -270.21392
2 -192.53361 -139.55831
3 -220.18590 -192.53361
4 -135.69672 -220.18590
5 -358.14644 -135.69672
6 -323.99318 -358.14644
7 -317.20746 -323.99318
8 -517.47561 -317.20746
9 -434.27615 -517.47561
10 -590.05029 -434.27615
11 -623.32771 -590.05029
12 336.12926 -623.32771
13 414.49450 336.12926
14 315.60624 414.49450
15 107.15925 315.60624
16 28.97130 107.15925
17 -33.38292 28.97130
18 -207.45597 -33.38292
19 -227.96325 -207.45597
20 -60.06602 -227.96325
21 -77.66295 -60.06602
22 -297.76769 -77.66295
23 -131.05404 -297.76769
24 -167.00770 -131.05404
25 -548.78955 -167.00770
26 -401.88533 -548.78955
27 -328.97921 -401.88533
28 32.80341 -328.97921
29 238.50366 32.80341
30 324.42271 238.50366
31 339.81203 324.42271
32 375.19404 339.81203
33 511.05121 375.19404
34 576.74494 511.05121
35 605.36438 576.74494
36 647.28621 605.36438
37 658.95390 647.28621
38 605.45336 658.95390
39 979.59718 605.45336
40 984.20407 979.59718
41 979.56395 984.20407
42 908.91016 979.56395
43 763.93361 908.91016
44 921.81365 763.93361
45 874.54327 921.81365
46 600.06280 874.54327
47 573.32232 600.06280
48 222.08904 573.32232
49 252.10404 222.08904
50 274.89992 252.10404
51 146.59483 274.89992
52 -100.77706 146.59483
53 153.31419 -100.77706
54 202.85523 153.31419
55 260.45952 202.85523
56 144.65070 260.45952
57 33.67203 144.65070
58 -288.98976 33.67203
59 -424.30496 -288.98976
60 -768.28289 -424.30496
61 -637.20458 -768.28289
62 -601.54059 -637.20458
63 -684.18615 -601.54059
64 -809.50500 -684.18615
65 -979.85244 -809.50500
66 -904.73895 -979.85244
67 -819.03445 -904.73895
68 -864.11676 -819.03445
69 -907.32740 -864.11676
> 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/7mfho1291107212.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/8mfho1291107212.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/9mfho1291107212.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/10x6y91291107212.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/110pef1291107212.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/1237v31291107212.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/13n2001291107212.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/14lirz1291107212.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/15p0qn1291107212.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/163ane1291107212.tab")
+ }
>
> try(system("convert tmp/185jx1291107212.ps tmp/185jx1291107212.png",intern=TRUE))
character(0)
> try(system("convert tmp/285jx1291107212.ps tmp/285jx1291107212.png",intern=TRUE))
character(0)
> try(system("convert tmp/3jei01291107212.ps tmp/3jei01291107212.png",intern=TRUE))
character(0)
> try(system("convert tmp/4jei01291107212.ps tmp/4jei01291107212.png",intern=TRUE))
character(0)
> try(system("convert tmp/5jei01291107212.ps tmp/5jei01291107212.png",intern=TRUE))
character(0)
> try(system("convert tmp/6u6il1291107212.ps tmp/6u6il1291107212.png",intern=TRUE))
character(0)
> try(system("convert tmp/7mfho1291107212.ps tmp/7mfho1291107212.png",intern=TRUE))
character(0)
> try(system("convert tmp/8mfho1291107212.ps tmp/8mfho1291107212.png",intern=TRUE))
character(0)
> try(system("convert tmp/9mfho1291107212.ps tmp/9mfho1291107212.png",intern=TRUE))
character(0)
> try(system("convert tmp/10x6y91291107212.ps tmp/10x6y91291107212.png",intern=TRUE))
character(0)
>
>
> proc.time()
user system elapsed
3.986 2.526 4.451