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(2981.85
+ ,10407
+ ,0.762253
+ ,14448.9
+ ,13953.3
+ ,3080.58
+ ,10463
+ ,0.768403
+ ,15023.9
+ ,14657.7
+ ,3106.22
+ ,10556
+ ,0.757518
+ ,17319.2
+ ,16686.2
+ ,3119.31
+ ,10646
+ ,0.772917
+ ,16080.7
+ ,15232.4
+ ,3061.26
+ ,10702
+ ,0.787774
+ ,15486.3
+ ,15014.1
+ ,3097.31
+ ,11353
+ ,0.82203
+ ,17046.4
+ ,16688.6
+ ,3161.69
+ ,11346
+ ,0.830772
+ ,14793.9
+ ,13969.6
+ ,3257.16
+ ,11451
+ ,0.813537
+ ,13666.7
+ ,14546.8
+ ,3277.01
+ ,11964
+ ,0.815927
+ ,17358.8
+ ,16292
+ ,3295.32
+ ,12574
+ ,0.832293
+ ,16091.8
+ ,15039
+ ,3363.99
+ ,13031
+ ,0.848464
+ ,17401.7
+ ,17433.8
+ ,3494.17
+ ,13812
+ ,0.843455
+ ,16467
+ ,17798.4
+ ,3667.03
+ ,14544
+ ,0.826241
+ ,16103.8
+ ,16870.9
+ ,3813.06
+ ,14931
+ ,0.837661
+ ,16422.6
+ ,16659.3
+ ,3917.96
+ ,14886
+ ,0.831947
+ ,19435.5
+ ,19620.4
+ ,3895.51
+ ,16005
+ ,0.81493
+ ,15810.1
+ ,15953.5
+ ,3801.06
+ ,17064
+ ,0.783085
+ ,17914.8
+ ,17420.9
+ ,3570.12
+ ,15168
+ ,0.790514
+ ,18197.2
+ ,17647.5
+ ,3701.61
+ ,16050
+ ,0.788395
+ ,16183.5
+ ,15200.8
+ ,3862.27
+ ,15839
+ ,0.780579
+ ,14781
+ ,15637.3
+ ,3970.1
+ ,15137
+ ,0.785731
+ ,18091.5
+ ,17124.5
+ ,4138.52
+ ,14954
+ ,0.792959
+ ,18318.8
+ ,17659.4
+ ,4199.75
+ ,15648
+ ,0.776337
+ ,18392.2
+ ,17815
+ ,4290.89
+ ,15305
+ ,0.75683
+ ,15952.5
+ ,16165.6
+ ,4443.91
+ ,15579
+ ,0.76929
+ ,17434.3
+ ,17416.6
+ ,4502.64
+ ,16348
+ ,0.764877
+ ,17214
+ ,16823.9
+ ,4356.98
+ ,15928
+ ,0.755173
+ ,19680.5
+ ,19171.2
+ ,4591.27
+ ,16171
+ ,0.739864
+ ,17216.8
+ ,16806.8
+ ,4696.96
+ ,15937
+ ,0.740138
+ ,18325.3
+ ,18112.8
+ ,4621.4
+ ,15713
+ ,0.745212
+ ,19303.5
+ ,18485.5
+ ,4562.84
+ ,15594
+ ,0.729076
+ ,18090.7
+ ,17668
+ ,4202.52
+ ,15683
+ ,0.734107
+ ,16166.3
+ ,16324.3
+ ,4296.49
+ ,16438
+ ,0.719632
+ ,18304.7
+ ,17877.5
+ ,4435.23
+ ,17032
+ ,0.702889
+ ,20380.1
+ ,20136.7
+ ,4105.18
+ ,17696
+ ,0.681013
+ ,18887.7
+ ,19307
+ ,4116.68
+ ,17745
+ ,0.686342
+ ,16316.5
+ ,17776.3
+ ,3844.49
+ ,19394
+ ,0.67944
+ ,18471.5
+ ,19861.3
+ ,3720.98
+ ,20148
+ ,0.678058
+ ,18754.9
+ ,18757
+ ,3674.4
+ ,20108
+ ,0.644039
+ ,18940.7
+ ,19879.3
+ ,3857.62
+ ,18584
+ ,0.63488
+ ,20228.5
+ ,21068.4
+ ,3801.06
+ ,18441
+ ,0.642797
+ ,19060.4
+ ,19358
+ ,3504.37
+ ,18391
+ ,0.642963
+ ,20262.9
+ ,20639.2
+ ,3032.6
+ ,19178
+ ,0.634115
+ ,19928.7
+ ,20008.1
+ ,3047.03
+ ,18079
+ ,0.66778
+ ,16058.8
+ ,18150.1
+ ,2962.34
+ ,18483
+ ,0.695894
+ ,20157.4
+ ,21180.4
+ ,2197.82
+ ,19644
+ ,0.750638
+ ,19663.3
+ ,20428.9
+ ,2014.45
+ ,19195
+ ,0.785423
+ ,15648.9
+ ,17241.2
+ ,1862.83
+ ,19650
+ ,0.74355
+ ,14380.5
+ ,15969.3
+ ,1905.41
+ ,20830
+ ,0.755344
+ ,13654.4
+ ,14972.4
+ ,1810.99
+ ,23595
+ ,0.782167
+ ,14085.9
+ ,14488.3
+ ,1670.07
+ ,22937
+ ,0.766284
+ ,15070.6
+ ,15885.1
+ ,1864.44
+ ,21814
+ ,0.75815
+ ,14206.9
+ ,14305.3
+ ,2052.02
+ ,21928
+ ,0.732601
+ ,13585.6
+ ,13891.5
+ ,2029.6
+ ,21777
+ ,0.71347
+ ,15413.2
+ ,15431.6
+ ,2070.83
+ ,21383
+ ,0.709824
+ ,14809.6
+ ,14199.3
+ ,2293.41
+ ,21467
+ ,0.700869
+ ,12625.3
+ ,13542.6
+ ,2443.27
+ ,22052
+ ,0.686719
+ ,16314.7
+ ,16226.3
+ ,2513.17
+ ,22680
+ ,0.674946
+ ,16045.9
+ ,16786.1
+ ,2466.92
+ ,24320
+ ,0.670511
+ ,16063.6
+ ,16034.3
+ ,2502.66
+ ,24977
+ ,0.684275
+ ,15851.3
+ ,16744.5
+ ,2539.91
+ ,25204
+ ,0.700673
+ ,14925.2
+ ,15955.4)
+ ,dim=c(5
+ ,61)
+ ,dimnames=list(c('BEL20'
+ ,'GoudkoersTeBrussel'
+ ,'EurosPerUSdollar'
+ ,'Uitvoer'
+ ,'Invoer')
+ ,1:61))
> y <- array(NA,dim=c(5,61),dimnames=list(c('BEL20','GoudkoersTeBrussel','EurosPerUSdollar','Uitvoer','Invoer'),1:61))
> 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
BEL20 GoudkoersTeBrussel EurosPerUSdollar Uitvoer Invoer
1 2981.85 10407 0.762253 14448.9 13953.3
2 3080.58 10463 0.768403 15023.9 14657.7
3 3106.22 10556 0.757518 17319.2 16686.2
4 3119.31 10646 0.772917 16080.7 15232.4
5 3061.26 10702 0.787774 15486.3 15014.1
6 3097.31 11353 0.822030 17046.4 16688.6
7 3161.69 11346 0.830772 14793.9 13969.6
8 3257.16 11451 0.813537 13666.7 14546.8
9 3277.01 11964 0.815927 17358.8 16292.0
10 3295.32 12574 0.832293 16091.8 15039.0
11 3363.99 13031 0.848464 17401.7 17433.8
12 3494.17 13812 0.843455 16467.0 17798.4
13 3667.03 14544 0.826241 16103.8 16870.9
14 3813.06 14931 0.837661 16422.6 16659.3
15 3917.96 14886 0.831947 19435.5 19620.4
16 3895.51 16005 0.814930 15810.1 15953.5
17 3801.06 17064 0.783085 17914.8 17420.9
18 3570.12 15168 0.790514 18197.2 17647.5
19 3701.61 16050 0.788395 16183.5 15200.8
20 3862.27 15839 0.780579 14781.0 15637.3
21 3970.10 15137 0.785731 18091.5 17124.5
22 4138.52 14954 0.792959 18318.8 17659.4
23 4199.75 15648 0.776337 18392.2 17815.0
24 4290.89 15305 0.756830 15952.5 16165.6
25 4443.91 15579 0.769290 17434.3 17416.6
26 4502.64 16348 0.764877 17214.0 16823.9
27 4356.98 15928 0.755173 19680.5 19171.2
28 4591.27 16171 0.739864 17216.8 16806.8
29 4696.96 15937 0.740138 18325.3 18112.8
30 4621.40 15713 0.745212 19303.5 18485.5
31 4562.84 15594 0.729076 18090.7 17668.0
32 4202.52 15683 0.734107 16166.3 16324.3
33 4296.49 16438 0.719632 18304.7 17877.5
34 4435.23 17032 0.702889 20380.1 20136.7
35 4105.18 17696 0.681013 18887.7 19307.0
36 4116.68 17745 0.686342 16316.5 17776.3
37 3844.49 19394 0.679440 18471.5 19861.3
38 3720.98 20148 0.678058 18754.9 18757.0
39 3674.40 20108 0.644039 18940.7 19879.3
40 3857.62 18584 0.634880 20228.5 21068.4
41 3801.06 18441 0.642797 19060.4 19358.0
42 3504.37 18391 0.642963 20262.9 20639.2
43 3032.60 19178 0.634115 19928.7 20008.1
44 3047.03 18079 0.667780 16058.8 18150.1
45 2962.34 18483 0.695894 20157.4 21180.4
46 2197.82 19644 0.750638 19663.3 20428.9
47 2014.45 19195 0.785423 15648.9 17241.2
48 1862.83 19650 0.743550 14380.5 15969.3
49 1905.41 20830 0.755344 13654.4 14972.4
50 1810.99 23595 0.782167 14085.9 14488.3
51 1670.07 22937 0.766284 15070.6 15885.1
52 1864.44 21814 0.758150 14206.9 14305.3
53 2052.02 21928 0.732601 13585.6 13891.5
54 2029.60 21777 0.713470 15413.2 15431.6
55 2070.83 21383 0.709824 14809.6 14199.3
56 2293.41 21467 0.700869 12625.3 13542.6
57 2443.27 22052 0.686719 16314.7 16226.3
58 2513.17 22680 0.674946 16045.9 16786.1
59 2466.92 24320 0.670511 16063.6 16034.3
60 2502.66 24977 0.684275 15851.3 16744.5
61 2539.91 25204 0.700673 14925.2 15955.4
> k <- length(x[1,])
> df <- as.data.frame(x)
> (mylm <- lm(df))
Call:
lm(formula = df)
Coefficients:
(Intercept) GoudkoersTeBrussel EurosPerUSdollar Uitvoer
2.127e+03 -9.908e-02 -1.106e+03 2.830e-01
Invoer
-6.092e-02
> (mysum <- summary(mylm))
Call:
lm(formula = df)
Residuals:
Min 1Q Median 3Q Max
-1472.88 -478.45 10.99 373.00 1036.33
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 2.127e+03 2.243e+03 0.948 0.347234
GoudkoersTeBrussel -9.908e-02 2.850e-02 -3.477 0.000988 ***
EurosPerUSdollar -1.106e+03 1.920e+03 -0.576 0.566810
Uitvoer 2.830e-01 1.189e-01 2.381 0.020705 *
Invoer -6.092e-02 1.195e-01 -0.510 0.612332
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 609.7 on 56 degrees of freedom
Multiple R-squared: 0.5337, Adjusted R-squared: 0.5004
F-statistic: 16.03 on 4 and 56 DF, p-value: 8.406e-09
> 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.001501e-03 4.003002e-03 9.979985e-01
[2,] 3.697755e-04 7.395510e-04 9.996302e-01
[3,] 9.680284e-05 1.936057e-04 9.999032e-01
[4,] 1.727967e-05 3.455933e-05 9.999827e-01
[5,] 2.166223e-06 4.332447e-06 9.999978e-01
[6,] 2.647120e-07 5.294241e-07 9.999997e-01
[7,] 1.066343e-07 2.132686e-07 9.999999e-01
[8,] 2.012746e-07 4.025491e-07 9.999998e-01
[9,] 6.330470e-08 1.266094e-07 9.999999e-01
[10,] 6.336119e-07 1.267224e-06 9.999994e-01
[11,] 3.316503e-07 6.633006e-07 9.999997e-01
[12,] 7.287868e-08 1.457574e-07 9.999999e-01
[13,] 1.852432e-08 3.704864e-08 1.000000e+00
[14,] 1.535528e-07 3.071056e-07 9.999998e-01
[15,] 2.460681e-06 4.921362e-06 9.999975e-01
[16,] 3.505271e-06 7.010543e-06 9.999965e-01
[17,] 1.266946e-05 2.533893e-05 9.999873e-01
[18,] 3.603223e-05 7.206445e-05 9.999640e-01
[19,] 5.600967e-05 1.120193e-04 9.999440e-01
[20,] 2.359746e-05 4.719491e-05 9.999764e-01
[21,] 2.841924e-05 5.683849e-05 9.999716e-01
[22,] 4.433185e-05 8.866370e-05 9.999557e-01
[23,] 4.369172e-05 8.738343e-05 9.999563e-01
[24,] 3.608875e-05 7.217750e-05 9.999639e-01
[25,] 3.734792e-05 7.469583e-05 9.999627e-01
[26,] 1.114718e-04 2.229436e-04 9.998885e-01
[27,] 1.315594e-03 2.631188e-03 9.986844e-01
[28,] 2.403165e-02 4.806330e-02 9.759684e-01
[29,] 1.668231e-01 3.336462e-01 8.331769e-01
[30,] 5.685294e-01 8.629413e-01 4.314706e-01
[31,] 9.616860e-01 7.662793e-02 3.831396e-02
[32,] 9.823730e-01 3.525404e-02 1.762702e-02
[33,] 9.862430e-01 2.751396e-02 1.375698e-02
[34,] 9.991442e-01 1.711634e-03 8.558172e-04
[35,] 9.997956e-01 4.087748e-04 2.043874e-04
[36,] 9.998315e-01 3.369033e-04 1.684517e-04
[37,] 9.997469e-01 5.061821e-04 2.530910e-04
[38,] 9.999284e-01 1.432061e-04 7.160303e-05
[39,] 9.999594e-01 8.126845e-05 4.063423e-05
[40,] 9.999985e-01 3.052818e-06 1.526409e-06
[41,] 9.999943e-01 1.142070e-05 5.710352e-06
[42,] 9.999731e-01 5.386551e-05 2.693276e-05
[43,] 9.999106e-01 1.787663e-04 8.938317e-05
[44,] 9.998600e-01 2.800398e-04 1.400199e-04
[45,] 9.989574e-01 2.085293e-03 1.042646e-03
[46,] 9.927785e-01 1.444305e-02 7.221526e-03
> postscript(file="/var/wessaorg/rcomp/tmp/1agn81353251529.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/2a1r91353251529.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/39c6c1353251529.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/418qz1353251529.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/56pg61353251529.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 = 61
Frequency = 1
1 2 3 4 5 6
-509.842457 -518.594621 -1021.858172 -720.841417 -601.970073 -803.078619
7 8 9 10 11 12
-257.823845 183.185722 -682.172578 -303.044680 -396.064355 92.722441
13 14 15 16 17 18
365.363187 459.247809 -119.001664 753.328267 222.265610 -254.439220
19 20 21 22 23 24
382.993642 937.651350 135.234106 261.769822 362.079095 987.697812
25 26 27 28 29 30
838.457023 994.744389 241.624105 1036.333332 884.957168 538.654650
31 32 33 34 35 36
743.918643 860.799061 502.939098 232.228957 315.627588 972.370240
37 38 39 40 41 42
373.004903 175.186951 102.794070 -167.171514 -2.725082 -566.486011
43 44 45 46 47 48
-913.922691 10.990307 -978.017908 -1472.879822 -720.231625 -591.570615
49 50 51 52 53 54
-274.246757 -216.657218 -633.955291 -411.634136 -90.379509 -572.377839
55 56 57 58 59 60
-478.449957 320.776810 -367.796703 -138.513960 -77.986256 141.429762
61
433.358675
> postscript(file="/var/wessaorg/rcomp/tmp/6htti1353251529.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 = 61
Frequency = 1
lag(myerror, k = 1) myerror
0 -509.842457 NA
1 -518.594621 -509.842457
2 -1021.858172 -518.594621
3 -720.841417 -1021.858172
4 -601.970073 -720.841417
5 -803.078619 -601.970073
6 -257.823845 -803.078619
7 183.185722 -257.823845
8 -682.172578 183.185722
9 -303.044680 -682.172578
10 -396.064355 -303.044680
11 92.722441 -396.064355
12 365.363187 92.722441
13 459.247809 365.363187
14 -119.001664 459.247809
15 753.328267 -119.001664
16 222.265610 753.328267
17 -254.439220 222.265610
18 382.993642 -254.439220
19 937.651350 382.993642
20 135.234106 937.651350
21 261.769822 135.234106
22 362.079095 261.769822
23 987.697812 362.079095
24 838.457023 987.697812
25 994.744389 838.457023
26 241.624105 994.744389
27 1036.333332 241.624105
28 884.957168 1036.333332
29 538.654650 884.957168
30 743.918643 538.654650
31 860.799061 743.918643
32 502.939098 860.799061
33 232.228957 502.939098
34 315.627588 232.228957
35 972.370240 315.627588
36 373.004903 972.370240
37 175.186951 373.004903
38 102.794070 175.186951
39 -167.171514 102.794070
40 -2.725082 -167.171514
41 -566.486011 -2.725082
42 -913.922691 -566.486011
43 10.990307 -913.922691
44 -978.017908 10.990307
45 -1472.879822 -978.017908
46 -720.231625 -1472.879822
47 -591.570615 -720.231625
48 -274.246757 -591.570615
49 -216.657218 -274.246757
50 -633.955291 -216.657218
51 -411.634136 -633.955291
52 -90.379509 -411.634136
53 -572.377839 -90.379509
54 -478.449957 -572.377839
55 320.776810 -478.449957
56 -367.796703 320.776810
57 -138.513960 -367.796703
58 -77.986256 -138.513960
59 141.429762 -77.986256
60 433.358675 141.429762
61 NA 433.358675
> dum1 <- dum[2:length(myerror),]
> dum1
lag(myerror, k = 1) myerror
[1,] -518.594621 -509.842457
[2,] -1021.858172 -518.594621
[3,] -720.841417 -1021.858172
[4,] -601.970073 -720.841417
[5,] -803.078619 -601.970073
[6,] -257.823845 -803.078619
[7,] 183.185722 -257.823845
[8,] -682.172578 183.185722
[9,] -303.044680 -682.172578
[10,] -396.064355 -303.044680
[11,] 92.722441 -396.064355
[12,] 365.363187 92.722441
[13,] 459.247809 365.363187
[14,] -119.001664 459.247809
[15,] 753.328267 -119.001664
[16,] 222.265610 753.328267
[17,] -254.439220 222.265610
[18,] 382.993642 -254.439220
[19,] 937.651350 382.993642
[20,] 135.234106 937.651350
[21,] 261.769822 135.234106
[22,] 362.079095 261.769822
[23,] 987.697812 362.079095
[24,] 838.457023 987.697812
[25,] 994.744389 838.457023
[26,] 241.624105 994.744389
[27,] 1036.333332 241.624105
[28,] 884.957168 1036.333332
[29,] 538.654650 884.957168
[30,] 743.918643 538.654650
[31,] 860.799061 743.918643
[32,] 502.939098 860.799061
[33,] 232.228957 502.939098
[34,] 315.627588 232.228957
[35,] 972.370240 315.627588
[36,] 373.004903 972.370240
[37,] 175.186951 373.004903
[38,] 102.794070 175.186951
[39,] -167.171514 102.794070
[40,] -2.725082 -167.171514
[41,] -566.486011 -2.725082
[42,] -913.922691 -566.486011
[43,] 10.990307 -913.922691
[44,] -978.017908 10.990307
[45,] -1472.879822 -978.017908
[46,] -720.231625 -1472.879822
[47,] -591.570615 -720.231625
[48,] -274.246757 -591.570615
[49,] -216.657218 -274.246757
[50,] -633.955291 -216.657218
[51,] -411.634136 -633.955291
[52,] -90.379509 -411.634136
[53,] -572.377839 -90.379509
[54,] -478.449957 -572.377839
[55,] 320.776810 -478.449957
[56,] -367.796703 320.776810
[57,] -138.513960 -367.796703
[58,] -77.986256 -138.513960
[59,] 141.429762 -77.986256
[60,] 433.358675 141.429762
> z <- as.data.frame(dum1)
> z
lag(myerror, k = 1) myerror
1 -518.594621 -509.842457
2 -1021.858172 -518.594621
3 -720.841417 -1021.858172
4 -601.970073 -720.841417
5 -803.078619 -601.970073
6 -257.823845 -803.078619
7 183.185722 -257.823845
8 -682.172578 183.185722
9 -303.044680 -682.172578
10 -396.064355 -303.044680
11 92.722441 -396.064355
12 365.363187 92.722441
13 459.247809 365.363187
14 -119.001664 459.247809
15 753.328267 -119.001664
16 222.265610 753.328267
17 -254.439220 222.265610
18 382.993642 -254.439220
19 937.651350 382.993642
20 135.234106 937.651350
21 261.769822 135.234106
22 362.079095 261.769822
23 987.697812 362.079095
24 838.457023 987.697812
25 994.744389 838.457023
26 241.624105 994.744389
27 1036.333332 241.624105
28 884.957168 1036.333332
29 538.654650 884.957168
30 743.918643 538.654650
31 860.799061 743.918643
32 502.939098 860.799061
33 232.228957 502.939098
34 315.627588 232.228957
35 972.370240 315.627588
36 373.004903 972.370240
37 175.186951 373.004903
38 102.794070 175.186951
39 -167.171514 102.794070
40 -2.725082 -167.171514
41 -566.486011 -2.725082
42 -913.922691 -566.486011
43 10.990307 -913.922691
44 -978.017908 10.990307
45 -1472.879822 -978.017908
46 -720.231625 -1472.879822
47 -591.570615 -720.231625
48 -274.246757 -591.570615
49 -216.657218 -274.246757
50 -633.955291 -216.657218
51 -411.634136 -633.955291
52 -90.379509 -411.634136
53 -572.377839 -90.379509
54 -478.449957 -572.377839
55 320.776810 -478.449957
56 -367.796703 320.776810
57 -138.513960 -367.796703
58 -77.986256 -138.513960
59 141.429762 -77.986256
60 433.358675 141.429762
> 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/7ceae1353251529.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/897to1353251529.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/9fuwe1353251529.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/10jeza1353251529.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/11mpg41353251529.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/127ehl1353251529.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/138rng1353251529.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/14gtrz1353251529.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/15nuc31353251529.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/167r641353251529.tab")
+ }
>
> try(system("convert tmp/1agn81353251529.ps tmp/1agn81353251529.png",intern=TRUE))
character(0)
> try(system("convert tmp/2a1r91353251529.ps tmp/2a1r91353251529.png",intern=TRUE))
character(0)
> try(system("convert tmp/39c6c1353251529.ps tmp/39c6c1353251529.png",intern=TRUE))
character(0)
> try(system("convert tmp/418qz1353251529.ps tmp/418qz1353251529.png",intern=TRUE))
character(0)
> try(system("convert tmp/56pg61353251529.ps tmp/56pg61353251529.png",intern=TRUE))
character(0)
> try(system("convert tmp/6htti1353251529.ps tmp/6htti1353251529.png",intern=TRUE))
character(0)
> try(system("convert tmp/7ceae1353251529.ps tmp/7ceae1353251529.png",intern=TRUE))
character(0)
> try(system("convert tmp/897to1353251529.ps tmp/897to1353251529.png",intern=TRUE))
character(0)
> try(system("convert tmp/9fuwe1353251529.ps tmp/9fuwe1353251529.png",intern=TRUE))
character(0)
> try(system("convert tmp/10jeza1353251529.ps tmp/10jeza1353251529.png",intern=TRUE))
character(0)
>
>
> proc.time()
user system elapsed
6.102 0.838 6.954