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(130
+ ,280
+ ,70
+ ,68.402973
+ ,3
+ ,15
+ ,135
+ ,120
+ ,33.983679
+ ,1
+ ,260
+ ,320
+ ,70
+ ,59.425505
+ ,3
+ ,200
+ ,158
+ ,110
+ ,34.384843
+ ,1
+ ,125
+ ,30
+ ,110
+ ,33.174094
+ ,2
+ ,200
+ ,125
+ ,90
+ ,49.120253
+ ,1
+ ,210
+ ,190
+ ,90
+ ,53.313813
+ ,1
+ ,220
+ ,35
+ ,120
+ ,18.042851
+ ,2
+ ,290
+ ,105
+ ,110
+ ,50.765
+ ,1
+ ,210
+ ,45
+ ,120
+ ,19.823573
+ ,2
+ ,140
+ ,105
+ ,110
+ ,40.400208
+ ,3
+ ,180
+ ,55
+ ,110
+ ,22.736446
+ ,2
+ ,280
+ ,25
+ ,110
+ ,41.445019
+ ,1
+ ,290
+ ,35
+ ,100
+ ,45.863324
+ ,1
+ ,90
+ ,20
+ ,110
+ ,35.782791
+ ,2
+ ,180
+ ,65
+ ,110
+ ,22.396513
+ ,1
+ ,80
+ ,25
+ ,100
+ ,64.533816
+ ,1
+ ,220
+ ,30
+ ,110
+ ,46.895644
+ ,3
+ ,190
+ ,80
+ ,100
+ ,44.330856
+ ,3
+ ,125
+ ,30
+ ,110
+ ,32.207582
+ ,2
+ ,200
+ ,25
+ ,110
+ ,31.435973
+ ,1
+ ,240
+ ,190
+ ,120
+ ,41.015492
+ ,3
+ ,135
+ ,25
+ ,110
+ ,28.025765
+ ,1
+ ,45
+ ,40
+ ,100
+ ,35.252444
+ ,1
+ ,280
+ ,45
+ ,110
+ ,23.804043
+ ,2
+ ,140
+ ,85
+ ,100
+ ,52.076897
+ ,3
+ ,170
+ ,90
+ ,110
+ ,53.371007
+ ,1
+ ,220
+ ,45
+ ,120
+ ,21.871292
+ ,2
+ ,250
+ ,90
+ ,110
+ ,31.072217
+ ,1
+ ,170
+ ,60
+ ,110
+ ,36.523683
+ ,1
+ ,260
+ ,40
+ ,110
+ ,39.241114
+ ,1
+ ,150
+ ,95
+ ,100
+ ,45.328074
+ ,2
+ ,180
+ ,55
+ ,110
+ ,26.734515
+ ,1
+ ,0
+ ,95
+ ,100
+ ,54.850917
+ ,1
+ ,220
+ ,90
+ ,100
+ ,40.105965
+ ,1
+ ,190
+ ,40
+ ,120
+ ,29.924285
+ ,1
+ ,170
+ ,120
+ ,130
+ ,30.450843
+ ,3
+ ,0
+ ,15
+ ,50
+ ,60.756112
+ ,3
+ ,0
+ ,50
+ ,50
+ ,63.005645
+ ,3
+ ,135
+ ,110
+ ,100
+ ,49.511874
+ ,3
+ ,0
+ ,110
+ ,100
+ ,50.828392
+ ,1
+ ,210
+ ,240
+ ,120
+ ,39.259197
+ ,2
+ ,140
+ ,140
+ ,100
+ ,3.7034
+ ,1
+ ,0
+ ,110
+ ,90
+ ,55.333142
+ ,3
+ ,240
+ ,30
+ ,110
+ ,41.998933
+ ,1
+ ,290
+ ,35
+ ,110
+ ,40.560159
+ ,1
+ ,0
+ ,95
+ ,80
+ ,68.235885
+ ,1
+ ,0
+ ,140
+ ,90
+ ,74.472949
+ ,1
+ ,0
+ ,120
+ ,90
+ ,72.801787
+ ,1
+ ,70
+ ,40
+ ,110
+ ,31.230054
+ ,2
+ ,230
+ ,55
+ ,110
+ ,53.131324
+ ,1
+ ,15
+ ,90
+ ,90
+ ,59.363993
+ ,2
+ ,200
+ ,35
+ ,110
+ ,38.839746
+ ,3
+ ,190
+ ,230
+ ,140
+ ,28.592785
+ ,3
+ ,200
+ ,110
+ ,100
+ ,46.658844
+ ,3
+ ,250
+ ,60
+ ,110
+ ,39.106174
+ ,1
+ ,140
+ ,25
+ ,110
+ ,27.753301
+ ,2
+ ,230
+ ,115
+ ,100
+ ,49.787445
+ ,1
+ ,200
+ ,110
+ ,100
+ ,51.592193
+ ,1
+ ,200
+ ,60
+ ,110
+ ,36.187559
+ ,1)
+ ,dim=c(5
+ ,60)
+ ,dimnames=list(c('fat'
+ ,'sugars'
+ ,'calories'
+ ,'rating'
+ ,'shelf')
+ ,1:60))
> y <- array(NA,dim=c(5,60),dimnames=list(c('fat','sugars','calories','rating','shelf'),1:60))
> 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 = '4'
> #'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
rating fat sugars calories shelf
1 68.40297 130 280 70 3
2 33.98368 15 135 120 1
3 59.42551 260 320 70 3
4 34.38484 200 158 110 1
5 33.17409 125 30 110 2
6 49.12025 200 125 90 1
7 53.31381 210 190 90 1
8 18.04285 220 35 120 2
9 50.76500 290 105 110 1
10 19.82357 210 45 120 2
11 40.40021 140 105 110 3
12 22.73645 180 55 110 2
13 41.44502 280 25 110 1
14 45.86332 290 35 100 1
15 35.78279 90 20 110 2
16 22.39651 180 65 110 1
17 64.53382 80 25 100 1
18 46.89564 220 30 110 3
19 44.33086 190 80 100 3
20 32.20758 125 30 110 2
21 31.43597 200 25 110 1
22 41.01549 240 190 120 3
23 28.02576 135 25 110 1
24 35.25244 45 40 100 1
25 23.80404 280 45 110 2
26 52.07690 140 85 100 3
27 53.37101 170 90 110 1
28 21.87129 220 45 120 2
29 31.07222 250 90 110 1
30 36.52368 170 60 110 1
31 39.24111 260 40 110 1
32 45.32807 150 95 100 2
33 26.73451 180 55 110 1
34 54.85092 0 95 100 1
35 40.10596 220 90 100 1
36 29.92429 190 40 120 1
37 30.45084 170 120 130 3
38 60.75611 0 15 50 3
39 63.00565 0 50 50 3
40 49.51187 135 110 100 3
41 50.82839 0 110 100 1
42 39.25920 210 240 120 2
43 3.70340 140 140 100 1
44 55.33314 0 110 90 3
45 41.99893 240 30 110 1
46 40.56016 290 35 110 1
47 68.23588 0 95 80 1
48 74.47295 0 140 90 1
49 72.80179 0 120 90 1
50 31.23005 70 40 110 2
51 53.13132 230 55 110 1
52 59.36399 15 90 90 2
53 38.83975 200 35 110 3
54 28.59278 190 230 140 3
55 46.65884 200 110 100 3
56 39.10617 250 60 110 1
57 27.75330 140 25 110 2
58 49.78744 230 115 100 1
59 51.59219 200 110 100 1
60 36.18756 200 60 110 1
> k <- length(x[1,])
> df <- as.data.frame(x)
> (mylm <- lm(df))
Call:
lm(formula = df)
Coefficients:
(Intercept) fat sugars calories shelf
99.74652 -0.02593 0.04900 -0.53383 -1.42358
> (mysum <- summary(mylm))
Call:
lm(formula = df)
Residuals:
Min 1Q Median 3Q Max
-44.468 -7.188 0.300 6.457 20.443
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 99.74652 10.53156 9.471 3.77e-13 ***
fat -0.02593 0.01750 -1.481 0.144
sugars 0.04900 0.02199 2.228 0.030 *
calories -0.53383 0.09920 -5.381 1.57e-06 ***
shelf -1.42358 1.69068 -0.842 0.403
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 10.57 on 55 degrees of freedom
Multiple R-squared: 0.5107, Adjusted R-squared: 0.4752
F-statistic: 14.35 on 4 and 55 DF, p-value: 4.361e-08
> 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,] 0.007866795 0.015733591 0.99213320
[2,] 0.350323651 0.700647303 0.64967635
[3,] 0.236758379 0.473516758 0.76324162
[4,] 0.247480524 0.494961048 0.75251948
[5,] 0.245173478 0.490346955 0.75482652
[6,] 0.206822901 0.413645803 0.79317710
[7,] 0.141295387 0.282590773 0.85870461
[8,] 0.091636656 0.183273312 0.90836334
[9,] 0.161240643 0.322481285 0.83875936
[10,] 0.308153292 0.616306585 0.69184671
[11,] 0.416335267 0.832670534 0.58366473
[12,] 0.331254762 0.662509523 0.66874524
[13,] 0.271357134 0.542714268 0.72864287
[14,] 0.223036123 0.446072246 0.77696388
[15,] 0.272791291 0.545582583 0.72720871
[16,] 0.247980926 0.495961852 0.75201907
[17,] 0.230487536 0.460975072 0.76951246
[18,] 0.229115787 0.458231574 0.77088421
[19,] 0.211540219 0.423080438 0.78845978
[20,] 0.279259272 0.558518543 0.72074073
[21,] 0.241730190 0.483460379 0.75826981
[22,] 0.198015515 0.396031029 0.80198449
[23,] 0.148401404 0.296802808 0.85159860
[24,] 0.111953989 0.223907979 0.88804601
[25,] 0.078149561 0.156299122 0.92185044
[26,] 0.078749728 0.157499455 0.92125027
[27,] 0.061174167 0.122348333 0.93882583
[28,] 0.042279000 0.084558001 0.95772100
[29,] 0.028644420 0.057288840 0.97135558
[30,] 0.019253208 0.038506415 0.98074679
[31,] 0.016373660 0.032747321 0.98362634
[32,] 0.014075267 0.028150534 0.98592473
[33,] 0.008938407 0.017876814 0.99106159
[34,] 0.005144263 0.010288527 0.99485574
[35,] 0.002738408 0.005476817 0.99726159
[36,] 0.974600394 0.050799212 0.02539961
[37,] 0.955858615 0.088282770 0.04414139
[38,] 0.934821788 0.130356424 0.06517821
[39,] 0.900454642 0.199090716 0.09954536
[40,] 0.871685862 0.256628276 0.12831414
[41,] 0.846662690 0.306674620 0.15333731
[42,] 0.868130147 0.263739705 0.13186985
[43,] 0.795796225 0.408407550 0.20420377
[44,] 0.959822710 0.080354581 0.04017729
[45,] 0.963303479 0.073393042 0.03669652
> postscript(file="/var/wessaorg/rcomp/tmp/1gmml1321894045.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/2pjew1321894045.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/3vr3t1321894045.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/4eeow1321894045.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/5s9gq1321894045.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 = 60
Frequency = 1
1 2 3 4 5 6
-0.05537405 -6.50648310 -7.62257963 -7.77434431 -3.23355985 -2.09838141
7 8 9 10 11 12
-0.83076832 -10.80856508 13.53632913 -9.77713582 2.12978603 -13.47035341
13 14 15 16 17 18
7.87734297 6.72660053 -1.04224365 -15.72389871 20.44264796 14.37454797
19 20 21 22 23 24
3.24354157 -4.20007185 -4.20578932 6.51068256 -9.30119244 -10.48118473
25 26 27 28 29 30
-9.32011659 9.44826266 13.76625462 -7.47015603 -6.45844909 -1.61097352
31 32 33 34 35 36
4.41986847 1.04508814 -10.89586476 5.25543899 -3.54075972 -1.37350977
37 38 39 40 41 42
2.89970801 -8.76333105 -8.22890989 5.52852938 0.49786606 0.10286508
43 44 45 46 47 48
-44.46757079 2.51150049 7.14919785 6.76171180 7.96385444 17.33405092
49 50 51 52 53 54
16.64295283 -7.09356613 16.79724818 6.48772622 5.55511242 1.50809636
55 56 57 58 59 60
4.36069450 3.04560377 -8.02044569 5.17490118 6.44688279 -1.16931516
> postscript(file="/var/wessaorg/rcomp/tmp/6t81z1321894045.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 = 60
Frequency = 1
lag(myerror, k = 1) myerror
0 -0.05537405 NA
1 -6.50648310 -0.05537405
2 -7.62257963 -6.50648310
3 -7.77434431 -7.62257963
4 -3.23355985 -7.77434431
5 -2.09838141 -3.23355985
6 -0.83076832 -2.09838141
7 -10.80856508 -0.83076832
8 13.53632913 -10.80856508
9 -9.77713582 13.53632913
10 2.12978603 -9.77713582
11 -13.47035341 2.12978603
12 7.87734297 -13.47035341
13 6.72660053 7.87734297
14 -1.04224365 6.72660053
15 -15.72389871 -1.04224365
16 20.44264796 -15.72389871
17 14.37454797 20.44264796
18 3.24354157 14.37454797
19 -4.20007185 3.24354157
20 -4.20578932 -4.20007185
21 6.51068256 -4.20578932
22 -9.30119244 6.51068256
23 -10.48118473 -9.30119244
24 -9.32011659 -10.48118473
25 9.44826266 -9.32011659
26 13.76625462 9.44826266
27 -7.47015603 13.76625462
28 -6.45844909 -7.47015603
29 -1.61097352 -6.45844909
30 4.41986847 -1.61097352
31 1.04508814 4.41986847
32 -10.89586476 1.04508814
33 5.25543899 -10.89586476
34 -3.54075972 5.25543899
35 -1.37350977 -3.54075972
36 2.89970801 -1.37350977
37 -8.76333105 2.89970801
38 -8.22890989 -8.76333105
39 5.52852938 -8.22890989
40 0.49786606 5.52852938
41 0.10286508 0.49786606
42 -44.46757079 0.10286508
43 2.51150049 -44.46757079
44 7.14919785 2.51150049
45 6.76171180 7.14919785
46 7.96385444 6.76171180
47 17.33405092 7.96385444
48 16.64295283 17.33405092
49 -7.09356613 16.64295283
50 16.79724818 -7.09356613
51 6.48772622 16.79724818
52 5.55511242 6.48772622
53 1.50809636 5.55511242
54 4.36069450 1.50809636
55 3.04560377 4.36069450
56 -8.02044569 3.04560377
57 5.17490118 -8.02044569
58 6.44688279 5.17490118
59 -1.16931516 6.44688279
60 NA -1.16931516
> dum1 <- dum[2:length(myerror),]
> dum1
lag(myerror, k = 1) myerror
[1,] -6.5064831 -0.05537405
[2,] -7.6225796 -6.50648310
[3,] -7.7743443 -7.62257963
[4,] -3.2335599 -7.77434431
[5,] -2.0983814 -3.23355985
[6,] -0.8307683 -2.09838141
[7,] -10.8085651 -0.83076832
[8,] 13.5363291 -10.80856508
[9,] -9.7771358 13.53632913
[10,] 2.1297860 -9.77713582
[11,] -13.4703534 2.12978603
[12,] 7.8773430 -13.47035341
[13,] 6.7266005 7.87734297
[14,] -1.0422437 6.72660053
[15,] -15.7238987 -1.04224365
[16,] 20.4426480 -15.72389871
[17,] 14.3745480 20.44264796
[18,] 3.2435416 14.37454797
[19,] -4.2000719 3.24354157
[20,] -4.2057893 -4.20007185
[21,] 6.5106826 -4.20578932
[22,] -9.3011924 6.51068256
[23,] -10.4811847 -9.30119244
[24,] -9.3201166 -10.48118473
[25,] 9.4482627 -9.32011659
[26,] 13.7662546 9.44826266
[27,] -7.4701560 13.76625462
[28,] -6.4584491 -7.47015603
[29,] -1.6109735 -6.45844909
[30,] 4.4198685 -1.61097352
[31,] 1.0450881 4.41986847
[32,] -10.8958648 1.04508814
[33,] 5.2554390 -10.89586476
[34,] -3.5407597 5.25543899
[35,] -1.3735098 -3.54075972
[36,] 2.8997080 -1.37350977
[37,] -8.7633311 2.89970801
[38,] -8.2289099 -8.76333105
[39,] 5.5285294 -8.22890989
[40,] 0.4978661 5.52852938
[41,] 0.1028651 0.49786606
[42,] -44.4675708 0.10286508
[43,] 2.5115005 -44.46757079
[44,] 7.1491978 2.51150049
[45,] 6.7617118 7.14919785
[46,] 7.9638544 6.76171180
[47,] 17.3340509 7.96385444
[48,] 16.6429528 17.33405092
[49,] -7.0935661 16.64295283
[50,] 16.7972482 -7.09356613
[51,] 6.4877262 16.79724818
[52,] 5.5551124 6.48772622
[53,] 1.5080964 5.55511242
[54,] 4.3606945 1.50809636
[55,] 3.0456038 4.36069450
[56,] -8.0204457 3.04560377
[57,] 5.1749012 -8.02044569
[58,] 6.4468828 5.17490118
[59,] -1.1693152 6.44688279
> z <- as.data.frame(dum1)
> z
lag(myerror, k = 1) myerror
1 -6.5064831 -0.05537405
2 -7.6225796 -6.50648310
3 -7.7743443 -7.62257963
4 -3.2335599 -7.77434431
5 -2.0983814 -3.23355985
6 -0.8307683 -2.09838141
7 -10.8085651 -0.83076832
8 13.5363291 -10.80856508
9 -9.7771358 13.53632913
10 2.1297860 -9.77713582
11 -13.4703534 2.12978603
12 7.8773430 -13.47035341
13 6.7266005 7.87734297
14 -1.0422437 6.72660053
15 -15.7238987 -1.04224365
16 20.4426480 -15.72389871
17 14.3745480 20.44264796
18 3.2435416 14.37454797
19 -4.2000719 3.24354157
20 -4.2057893 -4.20007185
21 6.5106826 -4.20578932
22 -9.3011924 6.51068256
23 -10.4811847 -9.30119244
24 -9.3201166 -10.48118473
25 9.4482627 -9.32011659
26 13.7662546 9.44826266
27 -7.4701560 13.76625462
28 -6.4584491 -7.47015603
29 -1.6109735 -6.45844909
30 4.4198685 -1.61097352
31 1.0450881 4.41986847
32 -10.8958648 1.04508814
33 5.2554390 -10.89586476
34 -3.5407597 5.25543899
35 -1.3735098 -3.54075972
36 2.8997080 -1.37350977
37 -8.7633311 2.89970801
38 -8.2289099 -8.76333105
39 5.5285294 -8.22890989
40 0.4978661 5.52852938
41 0.1028651 0.49786606
42 -44.4675708 0.10286508
43 2.5115005 -44.46757079
44 7.1491978 2.51150049
45 6.7617118 7.14919785
46 7.9638544 6.76171180
47 17.3340509 7.96385444
48 16.6429528 17.33405092
49 -7.0935661 16.64295283
50 16.7972482 -7.09356613
51 6.4877262 16.79724818
52 5.5551124 6.48772622
53 1.5080964 5.55511242
54 4.3606945 1.50809636
55 3.0456038 4.36069450
56 -8.0204457 3.04560377
57 5.1749012 -8.02044569
58 6.4468828 5.17490118
59 -1.1693152 6.44688279
> 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/73e8f1321894045.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/8pwn61321894045.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/9iar81321894045.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/10t65p1321894045.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/11x8jc1321894045.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/12m79f1321894045.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/13gp5a1321894045.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/14x7a01321894045.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/15liga1321894045.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/16f4gx1321894045.tab")
+ }
>
> try(system("convert tmp/1gmml1321894045.ps tmp/1gmml1321894045.png",intern=TRUE))
character(0)
> try(system("convert tmp/2pjew1321894045.ps tmp/2pjew1321894045.png",intern=TRUE))
character(0)
> try(system("convert tmp/3vr3t1321894045.ps tmp/3vr3t1321894045.png",intern=TRUE))
character(0)
> try(system("convert tmp/4eeow1321894045.ps tmp/4eeow1321894045.png",intern=TRUE))
character(0)
> try(system("convert tmp/5s9gq1321894045.ps tmp/5s9gq1321894045.png",intern=TRUE))
character(0)
> try(system("convert tmp/6t81z1321894045.ps tmp/6t81z1321894045.png",intern=TRUE))
character(0)
> try(system("convert tmp/73e8f1321894045.ps tmp/73e8f1321894045.png",intern=TRUE))
character(0)
> try(system("convert tmp/8pwn61321894045.ps tmp/8pwn61321894045.png",intern=TRUE))
character(0)
> try(system("convert tmp/9iar81321894045.ps tmp/9iar81321894045.png",intern=TRUE))
character(0)
> try(system("convert tmp/10t65p1321894045.ps tmp/10t65p1321894045.png",intern=TRUE))
character(0)
>
>
> proc.time()
user system elapsed
3.127 0.492 3.641