R version 2.12.0 (2010-10-15)
Copyright (C) 2010 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(104.37
+ ,1
+ ,167.16
+ ,101.56
+ ,100.93
+ ,104.89
+ ,2
+ ,179.84
+ ,102.13
+ ,101.18
+ ,105.15
+ ,3
+ ,174.44
+ ,102.39
+ ,101.11
+ ,105.72
+ ,4
+ ,180.35
+ ,102.42
+ ,102.42
+ ,106.38
+ ,5
+ ,193.17
+ ,103.87
+ ,102.37
+ ,106.40
+ ,6
+ ,195.16
+ ,104.44
+ ,101.95
+ ,106.47
+ ,7
+ ,202.43
+ ,104.97
+ ,102.20
+ ,106.59
+ ,8
+ ,189.91
+ ,105.17
+ ,103.35
+ ,106.76
+ ,9
+ ,195.98
+ ,105.35
+ ,103.65
+ ,107.35
+ ,10
+ ,212.09
+ ,104.65
+ ,102.06
+ ,107.81
+ ,11
+ ,205.81
+ ,106.62
+ ,102.66
+ ,108.03
+ ,12
+ ,204.31
+ ,107.05
+ ,102.32
+ ,109.08
+ ,1
+ ,196.07
+ ,112.30
+ ,102.21
+ ,109.86
+ ,2
+ ,199.98
+ ,114.70
+ ,102.33
+ ,110.29
+ ,3
+ ,199.1
+ ,115.40
+ ,104.41
+ ,110.34
+ ,4
+ ,198.31
+ ,115.64
+ ,104.33
+ ,110.59
+ ,5
+ ,195.72
+ ,115.66
+ ,105.27
+ ,110.64
+ ,6
+ ,223.04
+ ,114.50
+ ,105.34
+ ,110.83
+ ,7
+ ,238.41
+ ,115.14
+ ,104.88
+ ,111.51
+ ,8
+ ,259.73
+ ,115.41
+ ,105.49
+ ,113.32
+ ,9
+ ,326.54
+ ,119.32
+ ,105.90
+ ,115.89
+ ,10
+ ,335.15
+ ,124.77
+ ,105.39
+ ,116.51
+ ,11
+ ,321.81
+ ,130.96
+ ,104.40
+ ,117.44
+ ,12
+ ,368.62
+ ,141.02
+ ,106.19
+ ,118.25
+ ,1
+ ,369.59
+ ,150.60
+ ,106.54
+ ,118.65
+ ,2
+ ,425
+ ,151.10
+ ,108.26
+ ,118.52
+ ,3
+ ,439.72
+ ,157.19
+ ,106.95
+ ,119.07
+ ,4
+ ,362.23
+ ,157.28
+ ,108.32
+ ,119.12
+ ,5
+ ,328.76
+ ,156.54
+ ,108.35
+ ,119.28
+ ,6
+ ,348.55
+ ,159.62
+ ,109.29
+ ,119.30
+ ,7
+ ,328.18
+ ,163.77
+ ,109.46
+ ,119.44
+ ,8
+ ,329.34
+ ,165.08
+ ,109.50
+ ,119.57
+ ,9
+ ,295.55
+ ,164.75
+ ,109.84
+ ,119.93
+ ,10
+ ,237.38
+ ,163.93
+ ,108.73
+ ,120.03
+ ,11
+ ,226.85
+ ,157.51
+ ,109.38
+ ,119.66
+ ,12
+ ,220.14
+ ,153.36
+ ,109.97
+ ,119.46
+ ,1
+ ,239.36
+ ,156.83
+ ,111.10
+ ,119.48
+ ,2
+ ,224.69
+ ,154.98
+ ,110.53
+ ,119.56
+ ,3
+ ,230.98
+ ,155.02
+ ,110.23
+ ,119.43
+ ,4
+ ,233.47
+ ,153.34
+ ,109.41
+ ,119.57
+ ,5
+ ,256.7
+ ,153.19
+ ,108.94
+ ,119.59
+ ,6
+ ,253.41
+ ,152.80
+ ,109.81
+ ,119.50
+ ,7
+ ,224.95
+ ,152.97
+ ,109.20
+ ,119.54
+ ,8
+ ,210.37
+ ,152.96
+ ,109.45
+ ,119.56
+ ,9
+ ,191.09
+ ,152.35
+ ,110.61
+ ,119.61
+ ,10
+ ,198.85
+ ,151.88
+ ,109.44
+ ,119.64
+ ,11
+ ,211.04
+ ,150.27
+ ,109.77
+ ,119.60
+ ,12
+ ,206.25
+ ,148.80
+ ,108.04
+ ,119.71
+ ,1
+ ,201.19
+ ,149.28
+ ,109.65
+ ,119.72
+ ,2
+ ,194.37
+ ,148.64
+ ,111.69
+ ,119.66
+ ,3
+ ,191.08
+ ,150.36
+ ,111.65
+ ,119.76
+ ,4
+ ,192.87
+ ,149.69
+ ,112.04
+ ,119.80
+ ,5
+ ,181.61
+ ,152.94
+ ,111.42
+ ,119.88
+ ,6
+ ,157.67
+ ,155.18
+ ,112.25
+ ,119.78
+ ,7
+ ,196.14
+ ,156.32
+ ,111.46
+ ,120.08
+ ,8
+ ,246.35
+ ,156.25
+ ,111.62
+ ,120.22
+ ,9
+ ,271.9
+ ,155.52
+ ,111.77)
+ ,dim=c(5
+ ,57)
+ ,dimnames=list(c('Brood'
+ ,'Maand'
+ ,'Tarwe'
+ ,'Meel'
+ ,'Water')
+ ,1:57))
> y <- array(NA,dim=c(5,57),dimnames=list(c('Brood','Maand','Tarwe','Meel','Water'),1:57))
> 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'
> #'GNU S' R Code compiled by R2WASP v. 1.0.44 ()
> #Author: Prof. Dr. P. Wessa
> #To cite this work: AUTHOR(S), (YEAR), YOUR SOFTWARE TITLE (vNUMBER) in Free Statistics Software (v$_version), Office for Research Development and Education, URL http://www.wessa.net/rwasp_YOURPAGE.wasp/
> #Source of accompanying publication: Office for Research, Development, and Education
> #Technical description: Write here your technical program description (don't use hard returns!)
> library(lattice)
> library(lmtest)
Loading required package: zoo
> 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
Brood Maand Tarwe Meel Water
1 104.37 1 167.16 101.56 100.93
2 104.89 2 179.84 102.13 101.18
3 105.15 3 174.44 102.39 101.11
4 105.72 4 180.35 102.42 102.42
5 106.38 5 193.17 103.87 102.37
6 106.40 6 195.16 104.44 101.95
7 106.47 7 202.43 104.97 102.20
8 106.59 8 189.91 105.17 103.35
9 106.76 9 195.98 105.35 103.65
10 107.35 10 212.09 104.65 102.06
11 107.81 11 205.81 106.62 102.66
12 108.03 12 204.31 107.05 102.32
13 109.08 1 196.07 112.30 102.21
14 109.86 2 199.98 114.70 102.33
15 110.29 3 199.10 115.40 104.41
16 110.34 4 198.31 115.64 104.33
17 110.59 5 195.72 115.66 105.27
18 110.64 6 223.04 114.50 105.34
19 110.83 7 238.41 115.14 104.88
20 111.51 8 259.73 115.41 105.49
21 113.32 9 326.54 119.32 105.90
22 115.89 10 335.15 124.77 105.39
23 116.51 11 321.81 130.96 104.40
24 117.44 12 368.62 141.02 106.19
25 118.25 1 369.59 150.60 106.54
26 118.65 2 425.00 151.10 108.26
27 118.52 3 439.72 157.19 106.95
28 119.07 4 362.23 157.28 108.32
29 119.12 5 328.76 156.54 108.35
30 119.28 6 348.55 159.62 109.29
31 119.30 7 328.18 163.77 109.46
32 119.44 8 329.34 165.08 109.50
33 119.57 9 295.55 164.75 109.84
34 119.93 10 237.38 163.93 108.73
35 120.03 11 226.85 157.51 109.38
36 119.66 12 220.14 153.36 109.97
37 119.46 1 239.36 156.83 111.10
38 119.48 2 224.69 154.98 110.53
39 119.56 3 230.98 155.02 110.23
40 119.43 4 233.47 153.34 109.41
41 119.57 5 256.70 153.19 108.94
42 119.59 6 253.41 152.80 109.81
43 119.50 7 224.95 152.97 109.20
44 119.54 8 210.37 152.96 109.45
45 119.56 9 191.09 152.35 110.61
46 119.61 10 198.85 151.88 109.44
47 119.64 11 211.04 150.27 109.77
48 119.60 12 206.25 148.80 108.04
49 119.71 1 201.19 149.28 109.65
50 119.72 2 194.37 148.64 111.69
51 119.66 3 191.08 150.36 111.65
52 119.76 4 192.87 149.69 112.04
53 119.80 5 181.61 152.94 111.42
54 119.88 6 157.67 155.18 112.25
55 119.78 7 196.14 156.32 111.46
56 120.08 8 246.35 156.25 111.62
57 120.22 9 271.90 155.52 111.77
> k <- length(x[1,])
> df <- as.data.frame(x)
> (mylm <- lm(df))
Call:
lm(formula = df)
Coefficients:
(Intercept) Maand Tarwe Meel Water
26.936815 0.091968 0.006915 0.145336 0.617259
> (mysum <- summary(mylm))
Call:
lm(formula = df)
Residuals:
Min 1Q Median 3Q Max
-2.09191 -0.59464 0.05151 0.62312 2.86115
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 26.936815 10.936337 2.463 0.0171 *
Maand 0.091968 0.041760 2.202 0.0321 *
Tarwe 0.006915 0.002812 2.459 0.0173 *
Meel 0.145336 0.020765 6.999 4.97e-09 ***
Water 0.617259 0.123112 5.014 6.56e-06 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 1.062 on 52 degrees of freedom
Multiple R-squared: 0.967, Adjusted R-squared: 0.9645
F-statistic: 381.2 on 4 and 52 DF, p-value: < 2.2e-16
> if (n > n25) {
+ kp3 <- k + 3
+ nmkm3 <- n - k - 3
+ gqarr <- array(NA, dim=c(nmkm3-kp3+1,3))
+ numgqtests <- 0
+ numsignificant1 <- 0
+ numsignificant5 <- 0
+ numsignificant10 <- 0
+ for (mypoint in kp3:nmkm3) {
+ j <- 0
+ numgqtests <- numgqtests + 1
+ for (myalt in c('greater', 'two.sided', 'less')) {
+ j <- j + 1
+ gqarr[mypoint-kp3+1,j] <- gqtest(mylm, point=mypoint, alternative=myalt)$p.value
+ }
+ if (gqarr[mypoint-kp3+1,2] < 0.01) numsignificant1 <- numsignificant1 + 1
+ if (gqarr[mypoint-kp3+1,2] < 0.05) numsignificant5 <- numsignificant5 + 1
+ if (gqarr[mypoint-kp3+1,2] < 0.10) numsignificant10 <- numsignificant10 + 1
+ }
+ gqarr
+ }
[,1] [,2] [,3]
[1,] 1.795104e-02 3.590207e-02 9.820490e-01
[2,] 7.862910e-03 1.572582e-02 9.921371e-01
[3,] 1.822856e-03 3.645712e-03 9.981771e-01
[4,] 1.324154e-03 2.648309e-03 9.986758e-01
[5,] 5.263446e-04 1.052689e-03 9.994737e-01
[6,] 4.663652e-04 9.327303e-04 9.995336e-01
[7,] 1.290060e-04 2.580120e-04 9.998710e-01
[8,] 3.630009e-05 7.260017e-05 9.999637e-01
[9,] 1.099999e-05 2.199998e-05 9.999890e-01
[10,] 6.100619e-06 1.220124e-05 9.999939e-01
[11,] 4.338933e-06 8.677867e-06 9.999957e-01
[12,] 2.449998e-05 4.899996e-05 9.999755e-01
[13,] 8.388973e-04 1.677795e-03 9.991611e-01
[14,] 1.210034e-01 2.420067e-01 8.789966e-01
[15,] 1.759255e-01 3.518510e-01 8.240745e-01
[16,] 9.041563e-01 1.916874e-01 9.584370e-02
[17,] 1.000000e+00 5.684449e-10 2.842225e-10
[18,] 1.000000e+00 2.679290e-12 1.339645e-12
[19,] 1.000000e+00 7.142451e-13 3.571226e-13
[20,] 1.000000e+00 1.381885e-13 6.909423e-14
[21,] 1.000000e+00 2.788819e-13 1.394409e-13
[22,] 1.000000e+00 1.106892e-12 5.534458e-13
[23,] 1.000000e+00 2.459617e-12 1.229809e-12
[24,] 1.000000e+00 2.827602e-12 1.413801e-12
[25,] 1.000000e+00 4.078728e-12 2.039364e-12
[26,] 1.000000e+00 3.758618e-12 1.879309e-12
[27,] 1.000000e+00 6.492785e-12 3.246393e-12
[28,] 1.000000e+00 9.893608e-13 4.946804e-13
[29,] 1.000000e+00 4.199039e-12 2.099520e-12
[30,] 1.000000e+00 1.105223e-11 5.526114e-12
[31,] 1.000000e+00 4.551014e-11 2.275507e-11
[32,] 1.000000e+00 2.995639e-10 1.497820e-10
[33,] 1.000000e+00 1.068213e-09 5.341064e-10
[34,] 1.000000e+00 7.126603e-09 3.563301e-09
[35,] 1.000000e+00 1.891175e-08 9.455874e-09
[36,] 1.000000e+00 3.982058e-08 1.991029e-08
[37,] 1.000000e+00 9.256056e-08 4.628028e-08
[38,] 9.999998e-01 3.233542e-07 1.616771e-07
[39,] 9.999985e-01 2.938887e-06 1.469443e-06
[40,] 9.999869e-01 2.623880e-05 1.311940e-05
[41,] 9.998741e-01 2.518082e-04 1.259041e-04
[42,] 9.994298e-01 1.140402e-03 5.702008e-04
> postscript(file="/var/www/rcomp/tmp/1ts231292067142.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/rcomp/tmp/2412o1292067142.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/rcomp/tmp/3412o1292067142.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/rcomp/tmp/4412o1292067142.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/rcomp/tmp/5412o1292067142.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 = 57
Frequency = 1
1 2 3 4 5 6
-0.87497391 -0.77178151 -0.56098763 -0.93679268 -0.63728679 -0.54660904
7 8 9 10 11 12
-0.85019266 -1.47449877 -1.64977951 -0.17997286 -0.42518211 -0.13940436
13 14 15 16 17 18
1.28410810 1.52222378 0.48070760 0.45870232 0.05151441 -0.05399140
19 20 21 22 23 24
0.12867950 0.15351351 0.58820839 2.52942030 2.86115322 0.80851416
25 26 27 28 29 30
1.01509352 -0.19439161 -0.59463918 -0.45948306 -0.18097267 -1.27764906
31 32 33 34 35 36
-1.91683701 -2.09190759 -1.98212244 -0.50750812 0.10518066 -0.07142411
37 38 39 40 41 42
-0.59450143 0.05568413 0.17958431 0.69071473 0.89002148 0.36047012
43 44 45 46 47 48
0.72712543 0.62311773 0.05710735 0.75197887 0.63601220 1.81866900
49 50 51 52 53 54
1.91176114 0.71076152 0.35625579 0.20855416 0.14480691 -0.53949248
55 56 57
-0.67553254 -0.90329420 -1.01843565
> postscript(file="/var/www/rcomp/tmp/6xaj91292067142.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 = 57
Frequency = 1
lag(myerror, k = 1) myerror
0 -0.87497391 NA
1 -0.77178151 -0.87497391
2 -0.56098763 -0.77178151
3 -0.93679268 -0.56098763
4 -0.63728679 -0.93679268
5 -0.54660904 -0.63728679
6 -0.85019266 -0.54660904
7 -1.47449877 -0.85019266
8 -1.64977951 -1.47449877
9 -0.17997286 -1.64977951
10 -0.42518211 -0.17997286
11 -0.13940436 -0.42518211
12 1.28410810 -0.13940436
13 1.52222378 1.28410810
14 0.48070760 1.52222378
15 0.45870232 0.48070760
16 0.05151441 0.45870232
17 -0.05399140 0.05151441
18 0.12867950 -0.05399140
19 0.15351351 0.12867950
20 0.58820839 0.15351351
21 2.52942030 0.58820839
22 2.86115322 2.52942030
23 0.80851416 2.86115322
24 1.01509352 0.80851416
25 -0.19439161 1.01509352
26 -0.59463918 -0.19439161
27 -0.45948306 -0.59463918
28 -0.18097267 -0.45948306
29 -1.27764906 -0.18097267
30 -1.91683701 -1.27764906
31 -2.09190759 -1.91683701
32 -1.98212244 -2.09190759
33 -0.50750812 -1.98212244
34 0.10518066 -0.50750812
35 -0.07142411 0.10518066
36 -0.59450143 -0.07142411
37 0.05568413 -0.59450143
38 0.17958431 0.05568413
39 0.69071473 0.17958431
40 0.89002148 0.69071473
41 0.36047012 0.89002148
42 0.72712543 0.36047012
43 0.62311773 0.72712543
44 0.05710735 0.62311773
45 0.75197887 0.05710735
46 0.63601220 0.75197887
47 1.81866900 0.63601220
48 1.91176114 1.81866900
49 0.71076152 1.91176114
50 0.35625579 0.71076152
51 0.20855416 0.35625579
52 0.14480691 0.20855416
53 -0.53949248 0.14480691
54 -0.67553254 -0.53949248
55 -0.90329420 -0.67553254
56 -1.01843565 -0.90329420
57 NA -1.01843565
> dum1 <- dum[2:length(myerror),]
> dum1
lag(myerror, k = 1) myerror
[1,] -0.77178151 -0.87497391
[2,] -0.56098763 -0.77178151
[3,] -0.93679268 -0.56098763
[4,] -0.63728679 -0.93679268
[5,] -0.54660904 -0.63728679
[6,] -0.85019266 -0.54660904
[7,] -1.47449877 -0.85019266
[8,] -1.64977951 -1.47449877
[9,] -0.17997286 -1.64977951
[10,] -0.42518211 -0.17997286
[11,] -0.13940436 -0.42518211
[12,] 1.28410810 -0.13940436
[13,] 1.52222378 1.28410810
[14,] 0.48070760 1.52222378
[15,] 0.45870232 0.48070760
[16,] 0.05151441 0.45870232
[17,] -0.05399140 0.05151441
[18,] 0.12867950 -0.05399140
[19,] 0.15351351 0.12867950
[20,] 0.58820839 0.15351351
[21,] 2.52942030 0.58820839
[22,] 2.86115322 2.52942030
[23,] 0.80851416 2.86115322
[24,] 1.01509352 0.80851416
[25,] -0.19439161 1.01509352
[26,] -0.59463918 -0.19439161
[27,] -0.45948306 -0.59463918
[28,] -0.18097267 -0.45948306
[29,] -1.27764906 -0.18097267
[30,] -1.91683701 -1.27764906
[31,] -2.09190759 -1.91683701
[32,] -1.98212244 -2.09190759
[33,] -0.50750812 -1.98212244
[34,] 0.10518066 -0.50750812
[35,] -0.07142411 0.10518066
[36,] -0.59450143 -0.07142411
[37,] 0.05568413 -0.59450143
[38,] 0.17958431 0.05568413
[39,] 0.69071473 0.17958431
[40,] 0.89002148 0.69071473
[41,] 0.36047012 0.89002148
[42,] 0.72712543 0.36047012
[43,] 0.62311773 0.72712543
[44,] 0.05710735 0.62311773
[45,] 0.75197887 0.05710735
[46,] 0.63601220 0.75197887
[47,] 1.81866900 0.63601220
[48,] 1.91176114 1.81866900
[49,] 0.71076152 1.91176114
[50,] 0.35625579 0.71076152
[51,] 0.20855416 0.35625579
[52,] 0.14480691 0.20855416
[53,] -0.53949248 0.14480691
[54,] -0.67553254 -0.53949248
[55,] -0.90329420 -0.67553254
[56,] -1.01843565 -0.90329420
> z <- as.data.frame(dum1)
> z
lag(myerror, k = 1) myerror
1 -0.77178151 -0.87497391
2 -0.56098763 -0.77178151
3 -0.93679268 -0.56098763
4 -0.63728679 -0.93679268
5 -0.54660904 -0.63728679
6 -0.85019266 -0.54660904
7 -1.47449877 -0.85019266
8 -1.64977951 -1.47449877
9 -0.17997286 -1.64977951
10 -0.42518211 -0.17997286
11 -0.13940436 -0.42518211
12 1.28410810 -0.13940436
13 1.52222378 1.28410810
14 0.48070760 1.52222378
15 0.45870232 0.48070760
16 0.05151441 0.45870232
17 -0.05399140 0.05151441
18 0.12867950 -0.05399140
19 0.15351351 0.12867950
20 0.58820839 0.15351351
21 2.52942030 0.58820839
22 2.86115322 2.52942030
23 0.80851416 2.86115322
24 1.01509352 0.80851416
25 -0.19439161 1.01509352
26 -0.59463918 -0.19439161
27 -0.45948306 -0.59463918
28 -0.18097267 -0.45948306
29 -1.27764906 -0.18097267
30 -1.91683701 -1.27764906
31 -2.09190759 -1.91683701
32 -1.98212244 -2.09190759
33 -0.50750812 -1.98212244
34 0.10518066 -0.50750812
35 -0.07142411 0.10518066
36 -0.59450143 -0.07142411
37 0.05568413 -0.59450143
38 0.17958431 0.05568413
39 0.69071473 0.17958431
40 0.89002148 0.69071473
41 0.36047012 0.89002148
42 0.72712543 0.36047012
43 0.62311773 0.72712543
44 0.05710735 0.62311773
45 0.75197887 0.05710735
46 0.63601220 0.75197887
47 1.81866900 0.63601220
48 1.91176114 1.81866900
49 0.71076152 1.91176114
50 0.35625579 0.71076152
51 0.20855416 0.35625579
52 0.14480691 0.20855416
53 -0.53949248 0.14480691
54 -0.67553254 -0.53949248
55 -0.90329420 -0.67553254
56 -1.01843565 -0.90329420
> 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/rcomp/tmp/7pkiu1292067142.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/rcomp/tmp/8pkiu1292067142.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/rcomp/tmp/90bix1292067142.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/rcomp/tmp/100bix1292067142.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/rcomp/createtable file can be downloaded at http://www.wessa.net/cretab
> load(file="/var/www/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/rcomp/tmp/11mcg31292067142.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/rcomp/tmp/12pueq1292067142.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/rcomp/tmp/13wdbk1292067142.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/rcomp/tmp/146mbn1292067142.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/rcomp/tmp/152xu61292067143.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/rcomp/tmp/16g7ax1292067143.tab")
+ }
>
> try(system("convert tmp/1ts231292067142.ps tmp/1ts231292067142.png",intern=TRUE))
character(0)
> try(system("convert tmp/2412o1292067142.ps tmp/2412o1292067142.png",intern=TRUE))
character(0)
> try(system("convert tmp/3412o1292067142.ps tmp/3412o1292067142.png",intern=TRUE))
character(0)
> try(system("convert tmp/4412o1292067142.ps tmp/4412o1292067142.png",intern=TRUE))
character(0)
> try(system("convert tmp/5412o1292067142.ps tmp/5412o1292067142.png",intern=TRUE))
character(0)
> try(system("convert tmp/6xaj91292067142.ps tmp/6xaj91292067142.png",intern=TRUE))
character(0)
> try(system("convert tmp/7pkiu1292067142.ps tmp/7pkiu1292067142.png",intern=TRUE))
character(0)
> try(system("convert tmp/8pkiu1292067142.ps tmp/8pkiu1292067142.png",intern=TRUE))
character(0)
> try(system("convert tmp/90bix1292067142.ps tmp/90bix1292067142.png",intern=TRUE))
character(0)
> try(system("convert tmp/100bix1292067142.ps tmp/100bix1292067142.png",intern=TRUE))
character(0)
>
>
> proc.time()
user system elapsed
3.080 0.870 3.934