R version 2.13.0 (2011-04-13)
Copyright (C) 2011 The R Foundation for Statistical Computing
ISBN 3-900051-07-0
Platform: i486-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(68897,38683,44720,39525,45315,50380,40600,36279,42438,38064,31879,11379,70249,39253,47060,41697,38708,49267,39018,32228,40870,39383,34571,12066,70938,34077,45409,40809,37013,44953,37848,32745,39401,34931,33008,8620,68906,39556,50669,36432,40891,48428,36222,33425,39401,37967,34801,12657,69116,41519,51321,38529,41547,52073,38401,40898,40439,41888,37898,8771,68184,50530,47221,41756,45633,48138,39486,39341,41117,41629,29722,7054,56676,34870,35117,30169,30936,35699,33228,27733,33666,35429,27438,8170),dim=c(1,84),dimnames=list(c('Verkoop'),1:84))
> y <- array(NA,dim=c(1,84),dimnames=list(c('Verkoop'),1:84))
> for (i in 1:dim(x)[1])
+ {
+ for (j in 1:dim(x)[2])
+ {
+ y[i,j] <- as.numeric(x[i,j])
+ }
+ }
> par3 = 'Linear Trend'
> par2 = 'Include Monthly Dummies'
> par1 = '1'
> library(lattice)
> library(lmtest)
Loading required package: zoo
> 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
Verkoop M1 M2 M3 M4 M5 M6 M7 M8 M9 M10 M11 t
1 68897 1 0 0 0 0 0 0 0 0 0 0 1
2 38683 0 1 0 0 0 0 0 0 0 0 0 2
3 44720 0 0 1 0 0 0 0 0 0 0 0 3
4 39525 0 0 0 1 0 0 0 0 0 0 0 4
5 45315 0 0 0 0 1 0 0 0 0 0 0 5
6 50380 0 0 0 0 0 1 0 0 0 0 0 6
7 40600 0 0 0 0 0 0 1 0 0 0 0 7
8 36279 0 0 0 0 0 0 0 1 0 0 0 8
9 42438 0 0 0 0 0 0 0 0 1 0 0 9
10 38064 0 0 0 0 0 0 0 0 0 1 0 10
11 31879 0 0 0 0 0 0 0 0 0 0 1 11
12 11379 0 0 0 0 0 0 0 0 0 0 0 12
13 70249 1 0 0 0 0 0 0 0 0 0 0 13
14 39253 0 1 0 0 0 0 0 0 0 0 0 14
15 47060 0 0 1 0 0 0 0 0 0 0 0 15
16 41697 0 0 0 1 0 0 0 0 0 0 0 16
17 38708 0 0 0 0 1 0 0 0 0 0 0 17
18 49267 0 0 0 0 0 1 0 0 0 0 0 18
19 39018 0 0 0 0 0 0 1 0 0 0 0 19
20 32228 0 0 0 0 0 0 0 1 0 0 0 20
21 40870 0 0 0 0 0 0 0 0 1 0 0 21
22 39383 0 0 0 0 0 0 0 0 0 1 0 22
23 34571 0 0 0 0 0 0 0 0 0 0 1 23
24 12066 0 0 0 0 0 0 0 0 0 0 0 24
25 70938 1 0 0 0 0 0 0 0 0 0 0 25
26 34077 0 1 0 0 0 0 0 0 0 0 0 26
27 45409 0 0 1 0 0 0 0 0 0 0 0 27
28 40809 0 0 0 1 0 0 0 0 0 0 0 28
29 37013 0 0 0 0 1 0 0 0 0 0 0 29
30 44953 0 0 0 0 0 1 0 0 0 0 0 30
31 37848 0 0 0 0 0 0 1 0 0 0 0 31
32 32745 0 0 0 0 0 0 0 1 0 0 0 32
33 39401 0 0 0 0 0 0 0 0 1 0 0 33
34 34931 0 0 0 0 0 0 0 0 0 1 0 34
35 33008 0 0 0 0 0 0 0 0 0 0 1 35
36 8620 0 0 0 0 0 0 0 0 0 0 0 36
37 68906 1 0 0 0 0 0 0 0 0 0 0 37
38 39556 0 1 0 0 0 0 0 0 0 0 0 38
39 50669 0 0 1 0 0 0 0 0 0 0 0 39
40 36432 0 0 0 1 0 0 0 0 0 0 0 40
41 40891 0 0 0 0 1 0 0 0 0 0 0 41
42 48428 0 0 0 0 0 1 0 0 0 0 0 42
43 36222 0 0 0 0 0 0 1 0 0 0 0 43
44 33425 0 0 0 0 0 0 0 1 0 0 0 44
45 39401 0 0 0 0 0 0 0 0 1 0 0 45
46 37967 0 0 0 0 0 0 0 0 0 1 0 46
47 34801 0 0 0 0 0 0 0 0 0 0 1 47
48 12657 0 0 0 0 0 0 0 0 0 0 0 48
49 69116 1 0 0 0 0 0 0 0 0 0 0 49
50 41519 0 1 0 0 0 0 0 0 0 0 0 50
51 51321 0 0 1 0 0 0 0 0 0 0 0 51
52 38529 0 0 0 1 0 0 0 0 0 0 0 52
53 41547 0 0 0 0 1 0 0 0 0 0 0 53
54 52073 0 0 0 0 0 1 0 0 0 0 0 54
55 38401 0 0 0 0 0 0 1 0 0 0 0 55
56 40898 0 0 0 0 0 0 0 1 0 0 0 56
57 40439 0 0 0 0 0 0 0 0 1 0 0 57
58 41888 0 0 0 0 0 0 0 0 0 1 0 58
59 37898 0 0 0 0 0 0 0 0 0 0 1 59
60 8771 0 0 0 0 0 0 0 0 0 0 0 60
61 68184 1 0 0 0 0 0 0 0 0 0 0 61
62 50530 0 1 0 0 0 0 0 0 0 0 0 62
63 47221 0 0 1 0 0 0 0 0 0 0 0 63
64 41756 0 0 0 1 0 0 0 0 0 0 0 64
65 45633 0 0 0 0 1 0 0 0 0 0 0 65
66 48138 0 0 0 0 0 1 0 0 0 0 0 66
67 39486 0 0 0 0 0 0 1 0 0 0 0 67
68 39341 0 0 0 0 0 0 0 1 0 0 0 68
69 41117 0 0 0 0 0 0 0 0 1 0 0 69
70 41629 0 0 0 0 0 0 0 0 0 1 0 70
71 29722 0 0 0 0 0 0 0 0 0 0 1 71
72 7054 0 0 0 0 0 0 0 0 0 0 0 72
73 56676 1 0 0 0 0 0 0 0 0 0 0 73
74 34870 0 1 0 0 0 0 0 0 0 0 0 74
75 35117 0 0 1 0 0 0 0 0 0 0 0 75
76 30169 0 0 0 1 0 0 0 0 0 0 0 76
77 30936 0 0 0 0 1 0 0 0 0 0 0 77
78 35699 0 0 0 0 0 1 0 0 0 0 0 78
79 33228 0 0 0 0 0 0 1 0 0 0 0 79
80 27733 0 0 0 0 0 0 0 1 0 0 0 80
81 33666 0 0 0 0 0 0 0 0 1 0 0 81
82 35429 0 0 0 0 0 0 0 0 0 1 0 82
83 27438 0 0 0 0 0 0 0 0 0 0 1 83
84 8170 0 0 0 0 0 0 0 0 0 0 0 84
> k <- length(x[1,])
> df <- as.data.frame(x)
> (mylm <- lm(df))
Call:
lm(formula = df)
Coefficients:
(Intercept) M1 M2 M3 M4 M5
12478.07 57139.96 29412.84 35615.28 28156.44 29801.31
M6 M7 M8 M9 M10 M11
36841.76 27735.06 24625.65 29635.81 28542.54 22887.41
t
-55.44
> (mysum <- summary(mylm))
Call:
lm(formula = df)
Residuals:
Min 1Q Median 3Q Max
-9296.1 -2089.1 -209.8 2082.6 12076.7
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 12478.07 1741.67 7.164 5.86e-10 ***
M1 57139.96 2142.42 26.671 < 2e-16 ***
M2 29412.84 2140.80 13.739 < 2e-16 ***
M3 35615.28 2139.34 16.648 < 2e-16 ***
M4 28156.44 2138.03 13.169 < 2e-16 ***
M5 29801.31 2136.88 13.946 < 2e-16 ***
M6 36841.76 2135.88 17.249 < 2e-16 ***
M7 27735.06 2135.03 12.990 < 2e-16 ***
M8 24625.65 2134.34 11.538 < 2e-16 ***
M9 29635.81 2133.80 13.889 < 2e-16 ***
M10 28542.54 2133.41 13.379 < 2e-16 ***
M11 22887.41 2133.18 10.729 < 2e-16 ***
t -55.44 18.14 -3.056 0.00316 **
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 3991 on 71 degrees of freedom
Multiple R-squared: 0.9201, Adjusted R-squared: 0.9067
F-statistic: 68.18 on 12 and 71 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,] 3.146780e-03 6.293561e-03 0.9968532
[2,] 1.644854e-01 3.289707e-01 0.8355146
[3,] 8.254836e-02 1.650967e-01 0.9174516
[4,] 3.974874e-02 7.949747e-02 0.9602513
[5,] 2.997916e-02 5.995832e-02 0.9700208
[6,] 1.335339e-02 2.670677e-02 0.9866466
[7,] 7.014566e-03 1.402913e-02 0.9929854
[8,] 4.840564e-03 9.681128e-03 0.9951594
[9,] 2.077850e-03 4.155701e-03 0.9979221
[10,] 1.016845e-03 2.033690e-03 0.9989832
[11,] 2.216656e-03 4.433311e-03 0.9977833
[12,] 1.012165e-03 2.024330e-03 0.9989878
[13,] 4.452583e-04 8.905167e-04 0.9995547
[14,] 6.254207e-04 1.250841e-03 0.9993746
[15,] 6.087405e-04 1.217481e-03 0.9993913
[16,] 2.744310e-04 5.488620e-04 0.9997256
[17,] 1.514816e-04 3.029632e-04 0.9998485
[18,] 6.910386e-05 1.382077e-04 0.9999309
[19,] 6.963532e-05 1.392706e-04 0.9999304
[20,] 3.948923e-05 7.897847e-05 0.9999605
[21,] 2.530736e-05 5.061473e-05 0.9999747
[22,] 1.118279e-05 2.236559e-05 0.9999888
[23,] 3.017345e-05 6.034691e-05 0.9999698
[24,] 1.550072e-04 3.100145e-04 0.9998450
[25,] 1.487578e-04 2.975156e-04 0.9998512
[26,] 9.563115e-05 1.912623e-04 0.9999044
[27,] 5.128980e-05 1.025796e-04 0.9999487
[28,] 4.530621e-05 9.061241e-05 0.9999547
[29,] 6.205227e-05 1.241045e-04 0.9999379
[30,] 4.831522e-05 9.663044e-05 0.9999517
[31,] 1.079921e-04 2.159841e-04 0.9998920
[32,] 1.239747e-04 2.479495e-04 0.9998760
[33,] 1.273097e-04 2.546195e-04 0.9998727
[34,] 6.012454e-05 1.202491e-04 0.9999399
[35,] 2.784980e-04 5.569960e-04 0.9997215
[36,] 4.232211e-04 8.464421e-04 0.9995768
[37,] 3.497204e-04 6.994408e-04 0.9996503
[38,] 2.966902e-04 5.933805e-04 0.9997033
[39,] 2.913777e-04 5.827555e-04 0.9997086
[40,] 3.416710e-04 6.833421e-04 0.9996583
[41,] 9.410651e-04 1.882130e-03 0.9990589
[42,] 1.140517e-03 2.281033e-03 0.9988595
[43,] 2.711609e-03 5.423218e-03 0.9972884
[44,] 2.145868e-03 4.291736e-03 0.9978541
[45,] 3.311774e-02 6.623548e-02 0.9668823
[46,] 2.330985e-02 4.661970e-02 0.9766901
[47,] 1.371695e-01 2.743389e-01 0.8628305
[48,] 1.176277e-01 2.352555e-01 0.8823723
[49,] 9.716055e-02 1.943211e-01 0.9028395
[50,] 1.983265e-01 3.966529e-01 0.8016735
[51,] 2.938834e-01 5.877667e-01 0.7061166
[52,] 1.890496e-01 3.780992e-01 0.8109504
[53,] 3.938027e-01 7.876054e-01 0.6061973
> postscript(file="/var/wessaorg/rcomp/tmp/1glk21322262163.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/2v9y81322262163.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/31yuk1322262163.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/4xp9k1322262163.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/5ylew1322262163.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 = 84
Frequency = 1
1 2 3 4 5 6
-665.58929 -3097.01786 -3207.01786 -887.73214 3312.83929 1392.83929
7 8 9 10 11 12
774.98214 -381.16071 823.12500 -2402.16071 -2876.58929 -433.73214
13 14 15 16 17 18
1351.75000 -1861.67857 -201.67857 1949.60714 -2628.82143 945.17857
19 20 21 22 23 24
-141.67857 -3766.82143 -79.53571 -417.82143 480.75000 918.60714
25 26 27 28 29 30
2706.08929 -6372.33929 -1187.33929 1726.94643 -3658.48214 -2703.48214
31 32 33 34 35 36
-646.33929 -2584.48214 -883.19643 -4204.48214 -416.91071 -1862.05357
37 38 39 40 41 42
1339.42857 -228.00000 4738.00000 -1984.71429 884.85714 1436.85714
43 44 45 46 47 48
-1607.00000 -1239.14286 -217.85714 -503.14286 2041.42857 2840.28571
49 50 51 52 53 54
2214.76786 2400.33929 6055.33929 777.62500 2206.19643 5747.19643
55 56 57 58 59 60
1237.33929 6899.19643 1485.48214 4083.19643 5803.76786 -380.37500
61 62 63 64 65 66
1948.10714 12076.67857 2620.67857 4669.96429 6957.53571 2477.53571
67 68 69 70 71 72
2987.67857 6007.53571 2828.82143 4489.53571 -1706.89286 -1432.03571
73 74 75 76 77 78
-8894.55357 -2917.98214 -8817.98214 -6251.69643 -7074.12500 -9296.12500
79 80 81 82 83 84
-2604.98214 -4935.12500 -3956.83929 -1045.12500 -3325.55357 349.30357
> postscript(file="/var/wessaorg/rcomp/tmp/6kss61322262163.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 = 84
Frequency = 1
lag(myerror, k = 1) myerror
0 -665.58929 NA
1 -3097.01786 -665.58929
2 -3207.01786 -3097.01786
3 -887.73214 -3207.01786
4 3312.83929 -887.73214
5 1392.83929 3312.83929
6 774.98214 1392.83929
7 -381.16071 774.98214
8 823.12500 -381.16071
9 -2402.16071 823.12500
10 -2876.58929 -2402.16071
11 -433.73214 -2876.58929
12 1351.75000 -433.73214
13 -1861.67857 1351.75000
14 -201.67857 -1861.67857
15 1949.60714 -201.67857
16 -2628.82143 1949.60714
17 945.17857 -2628.82143
18 -141.67857 945.17857
19 -3766.82143 -141.67857
20 -79.53571 -3766.82143
21 -417.82143 -79.53571
22 480.75000 -417.82143
23 918.60714 480.75000
24 2706.08929 918.60714
25 -6372.33929 2706.08929
26 -1187.33929 -6372.33929
27 1726.94643 -1187.33929
28 -3658.48214 1726.94643
29 -2703.48214 -3658.48214
30 -646.33929 -2703.48214
31 -2584.48214 -646.33929
32 -883.19643 -2584.48214
33 -4204.48214 -883.19643
34 -416.91071 -4204.48214
35 -1862.05357 -416.91071
36 1339.42857 -1862.05357
37 -228.00000 1339.42857
38 4738.00000 -228.00000
39 -1984.71429 4738.00000
40 884.85714 -1984.71429
41 1436.85714 884.85714
42 -1607.00000 1436.85714
43 -1239.14286 -1607.00000
44 -217.85714 -1239.14286
45 -503.14286 -217.85714
46 2041.42857 -503.14286
47 2840.28571 2041.42857
48 2214.76786 2840.28571
49 2400.33929 2214.76786
50 6055.33929 2400.33929
51 777.62500 6055.33929
52 2206.19643 777.62500
53 5747.19643 2206.19643
54 1237.33929 5747.19643
55 6899.19643 1237.33929
56 1485.48214 6899.19643
57 4083.19643 1485.48214
58 5803.76786 4083.19643
59 -380.37500 5803.76786
60 1948.10714 -380.37500
61 12076.67857 1948.10714
62 2620.67857 12076.67857
63 4669.96429 2620.67857
64 6957.53571 4669.96429
65 2477.53571 6957.53571
66 2987.67857 2477.53571
67 6007.53571 2987.67857
68 2828.82143 6007.53571
69 4489.53571 2828.82143
70 -1706.89286 4489.53571
71 -1432.03571 -1706.89286
72 -8894.55357 -1432.03571
73 -2917.98214 -8894.55357
74 -8817.98214 -2917.98214
75 -6251.69643 -8817.98214
76 -7074.12500 -6251.69643
77 -9296.12500 -7074.12500
78 -2604.98214 -9296.12500
79 -4935.12500 -2604.98214
80 -3956.83929 -4935.12500
81 -1045.12500 -3956.83929
82 -3325.55357 -1045.12500
83 349.30357 -3325.55357
84 NA 349.30357
> dum1 <- dum[2:length(myerror),]
> dum1
lag(myerror, k = 1) myerror
[1,] -3097.01786 -665.58929
[2,] -3207.01786 -3097.01786
[3,] -887.73214 -3207.01786
[4,] 3312.83929 -887.73214
[5,] 1392.83929 3312.83929
[6,] 774.98214 1392.83929
[7,] -381.16071 774.98214
[8,] 823.12500 -381.16071
[9,] -2402.16071 823.12500
[10,] -2876.58929 -2402.16071
[11,] -433.73214 -2876.58929
[12,] 1351.75000 -433.73214
[13,] -1861.67857 1351.75000
[14,] -201.67857 -1861.67857
[15,] 1949.60714 -201.67857
[16,] -2628.82143 1949.60714
[17,] 945.17857 -2628.82143
[18,] -141.67857 945.17857
[19,] -3766.82143 -141.67857
[20,] -79.53571 -3766.82143
[21,] -417.82143 -79.53571
[22,] 480.75000 -417.82143
[23,] 918.60714 480.75000
[24,] 2706.08929 918.60714
[25,] -6372.33929 2706.08929
[26,] -1187.33929 -6372.33929
[27,] 1726.94643 -1187.33929
[28,] -3658.48214 1726.94643
[29,] -2703.48214 -3658.48214
[30,] -646.33929 -2703.48214
[31,] -2584.48214 -646.33929
[32,] -883.19643 -2584.48214
[33,] -4204.48214 -883.19643
[34,] -416.91071 -4204.48214
[35,] -1862.05357 -416.91071
[36,] 1339.42857 -1862.05357
[37,] -228.00000 1339.42857
[38,] 4738.00000 -228.00000
[39,] -1984.71429 4738.00000
[40,] 884.85714 -1984.71429
[41,] 1436.85714 884.85714
[42,] -1607.00000 1436.85714
[43,] -1239.14286 -1607.00000
[44,] -217.85714 -1239.14286
[45,] -503.14286 -217.85714
[46,] 2041.42857 -503.14286
[47,] 2840.28571 2041.42857
[48,] 2214.76786 2840.28571
[49,] 2400.33929 2214.76786
[50,] 6055.33929 2400.33929
[51,] 777.62500 6055.33929
[52,] 2206.19643 777.62500
[53,] 5747.19643 2206.19643
[54,] 1237.33929 5747.19643
[55,] 6899.19643 1237.33929
[56,] 1485.48214 6899.19643
[57,] 4083.19643 1485.48214
[58,] 5803.76786 4083.19643
[59,] -380.37500 5803.76786
[60,] 1948.10714 -380.37500
[61,] 12076.67857 1948.10714
[62,] 2620.67857 12076.67857
[63,] 4669.96429 2620.67857
[64,] 6957.53571 4669.96429
[65,] 2477.53571 6957.53571
[66,] 2987.67857 2477.53571
[67,] 6007.53571 2987.67857
[68,] 2828.82143 6007.53571
[69,] 4489.53571 2828.82143
[70,] -1706.89286 4489.53571
[71,] -1432.03571 -1706.89286
[72,] -8894.55357 -1432.03571
[73,] -2917.98214 -8894.55357
[74,] -8817.98214 -2917.98214
[75,] -6251.69643 -8817.98214
[76,] -7074.12500 -6251.69643
[77,] -9296.12500 -7074.12500
[78,] -2604.98214 -9296.12500
[79,] -4935.12500 -2604.98214
[80,] -3956.83929 -4935.12500
[81,] -1045.12500 -3956.83929
[82,] -3325.55357 -1045.12500
[83,] 349.30357 -3325.55357
> z <- as.data.frame(dum1)
> z
lag(myerror, k = 1) myerror
1 -3097.01786 -665.58929
2 -3207.01786 -3097.01786
3 -887.73214 -3207.01786
4 3312.83929 -887.73214
5 1392.83929 3312.83929
6 774.98214 1392.83929
7 -381.16071 774.98214
8 823.12500 -381.16071
9 -2402.16071 823.12500
10 -2876.58929 -2402.16071
11 -433.73214 -2876.58929
12 1351.75000 -433.73214
13 -1861.67857 1351.75000
14 -201.67857 -1861.67857
15 1949.60714 -201.67857
16 -2628.82143 1949.60714
17 945.17857 -2628.82143
18 -141.67857 945.17857
19 -3766.82143 -141.67857
20 -79.53571 -3766.82143
21 -417.82143 -79.53571
22 480.75000 -417.82143
23 918.60714 480.75000
24 2706.08929 918.60714
25 -6372.33929 2706.08929
26 -1187.33929 -6372.33929
27 1726.94643 -1187.33929
28 -3658.48214 1726.94643
29 -2703.48214 -3658.48214
30 -646.33929 -2703.48214
31 -2584.48214 -646.33929
32 -883.19643 -2584.48214
33 -4204.48214 -883.19643
34 -416.91071 -4204.48214
35 -1862.05357 -416.91071
36 1339.42857 -1862.05357
37 -228.00000 1339.42857
38 4738.00000 -228.00000
39 -1984.71429 4738.00000
40 884.85714 -1984.71429
41 1436.85714 884.85714
42 -1607.00000 1436.85714
43 -1239.14286 -1607.00000
44 -217.85714 -1239.14286
45 -503.14286 -217.85714
46 2041.42857 -503.14286
47 2840.28571 2041.42857
48 2214.76786 2840.28571
49 2400.33929 2214.76786
50 6055.33929 2400.33929
51 777.62500 6055.33929
52 2206.19643 777.62500
53 5747.19643 2206.19643
54 1237.33929 5747.19643
55 6899.19643 1237.33929
56 1485.48214 6899.19643
57 4083.19643 1485.48214
58 5803.76786 4083.19643
59 -380.37500 5803.76786
60 1948.10714 -380.37500
61 12076.67857 1948.10714
62 2620.67857 12076.67857
63 4669.96429 2620.67857
64 6957.53571 4669.96429
65 2477.53571 6957.53571
66 2987.67857 2477.53571
67 6007.53571 2987.67857
68 2828.82143 6007.53571
69 4489.53571 2828.82143
70 -1706.89286 4489.53571
71 -1432.03571 -1706.89286
72 -8894.55357 -1432.03571
73 -2917.98214 -8894.55357
74 -8817.98214 -2917.98214
75 -6251.69643 -8817.98214
76 -7074.12500 -6251.69643
77 -9296.12500 -7074.12500
78 -2604.98214 -9296.12500
79 -4935.12500 -2604.98214
80 -3956.83929 -4935.12500
81 -1045.12500 -3956.83929
82 -3325.55357 -1045.12500
83 349.30357 -3325.55357
> 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/782iv1322262163.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/8t24a1322262163.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/9nocj1322262163.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/10ybrz1322262163.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/119mj71322262163.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/12btsa1322262163.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/13f67z1322262163.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/14z0671322262164.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/152o3w1322262164.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/16q86e1322262164.tab")
+ }
>
> try(system("convert tmp/1glk21322262163.ps tmp/1glk21322262163.png",intern=TRUE))
character(0)
> try(system("convert tmp/2v9y81322262163.ps tmp/2v9y81322262163.png",intern=TRUE))
character(0)
> try(system("convert tmp/31yuk1322262163.ps tmp/31yuk1322262163.png",intern=TRUE))
character(0)
> try(system("convert tmp/4xp9k1322262163.ps tmp/4xp9k1322262163.png",intern=TRUE))
character(0)
> try(system("convert tmp/5ylew1322262163.ps tmp/5ylew1322262163.png",intern=TRUE))
character(0)
> try(system("convert tmp/6kss61322262163.ps tmp/6kss61322262163.png",intern=TRUE))
character(0)
> try(system("convert tmp/782iv1322262163.ps tmp/782iv1322262163.png",intern=TRUE))
character(0)
> try(system("convert tmp/8t24a1322262163.ps tmp/8t24a1322262163.png",intern=TRUE))
character(0)
> try(system("convert tmp/9nocj1322262163.ps tmp/9nocj1322262163.png",intern=TRUE))
character(0)
> try(system("convert tmp/10ybrz1322262163.ps tmp/10ybrz1322262163.png",intern=TRUE))
character(0)
>
>
> proc.time()
user system elapsed
3.483 0.615 4.222