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
+ ,15
+ ,135
+ ,120
+ ,33.983679
+ ,260
+ ,320
+ ,70
+ ,59.425505
+ ,200
+ ,158
+ ,110
+ ,34.384843
+ ,125
+ ,30
+ ,110
+ ,33.174094
+ ,200
+ ,125
+ ,90
+ ,49.120253
+ ,210
+ ,190
+ ,90
+ ,53.313813
+ ,220
+ ,35
+ ,120
+ ,18.042851
+ ,290
+ ,105
+ ,110
+ ,50.765
+ ,210
+ ,45
+ ,120
+ ,19.823573
+ ,140
+ ,105
+ ,110
+ ,40.400208
+ ,180
+ ,55
+ ,110
+ ,22.736446
+ ,280
+ ,25
+ ,110
+ ,41.445019
+ ,290
+ ,35
+ ,100
+ ,45.863324
+ ,90
+ ,20
+ ,110
+ ,35.782791
+ ,180
+ ,65
+ ,110
+ ,22.396513
+ ,80
+ ,25
+ ,100
+ ,64.533816
+ ,220
+ ,30
+ ,110
+ ,46.895644
+ ,190
+ ,80
+ ,100
+ ,44.330856
+ ,125
+ ,30
+ ,110
+ ,32.207582
+ ,200
+ ,25
+ ,110
+ ,31.435973
+ ,240
+ ,190
+ ,120
+ ,41.015492
+ ,135
+ ,25
+ ,110
+ ,28.025765
+ ,45
+ ,40
+ ,100
+ ,35.252444
+ ,280
+ ,45
+ ,110
+ ,23.804043
+ ,140
+ ,85
+ ,100
+ ,52.076897
+ ,170
+ ,90
+ ,110
+ ,53.371007
+ ,220
+ ,45
+ ,120
+ ,21.871292
+ ,250
+ ,90
+ ,110
+ ,31.072217
+ ,170
+ ,60
+ ,110
+ ,36.523683
+ ,260
+ ,40
+ ,110
+ ,39.241114
+ ,150
+ ,95
+ ,100
+ ,45.328074
+ ,180
+ ,55
+ ,110
+ ,26.734515
+ ,0
+ ,95
+ ,100
+ ,54.850917
+ ,220
+ ,90
+ ,100
+ ,40.105965
+ ,190
+ ,40
+ ,120
+ ,29.924285
+ ,170
+ ,120
+ ,130
+ ,30.450843
+ ,0
+ ,15
+ ,50
+ ,60.756112
+ ,0
+ ,50
+ ,50
+ ,63.005645
+ ,135
+ ,110
+ ,100
+ ,49.511874
+ ,0
+ ,110
+ ,100
+ ,50.828392
+ ,210
+ ,240
+ ,120
+ ,39.259197
+ ,140
+ ,140
+ ,100
+ ,3.7034
+ ,0
+ ,110
+ ,90
+ ,55.333142
+ ,240
+ ,30
+ ,110
+ ,41.998933
+ ,290
+ ,35
+ ,110
+ ,40.560159
+ ,0
+ ,95
+ ,80
+ ,68.235885
+ ,0
+ ,140
+ ,90
+ ,74.472949
+ ,0
+ ,120
+ ,90
+ ,72.801787
+ ,70
+ ,40
+ ,110
+ ,31.230054
+ ,230
+ ,55
+ ,110
+ ,53.131324
+ ,15
+ ,90
+ ,90
+ ,59.363993
+ ,200
+ ,35
+ ,110
+ ,38.839746
+ ,190
+ ,230
+ ,140
+ ,28.592785
+ ,200
+ ,110
+ ,100
+ ,46.658844
+ ,250
+ ,60
+ ,110
+ ,39.106174
+ ,140
+ ,25
+ ,110
+ ,27.753301
+ ,230
+ ,115
+ ,100
+ ,49.787445
+ ,200
+ ,110
+ ,100
+ ,51.592193
+ ,200
+ ,60
+ ,110
+ ,36.187559)
+ ,dim=c(4
+ ,60)
+ ,dimnames=list(c('fat'
+ ,'sugars'
+ ,'calories'
+ ,'rating')
+ ,1:60))
> y <- array(NA,dim=c(4,60),dimnames=list(c('fat','sugars','calories','rating'),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
1 68.40297 130 280 70
2 33.98368 15 135 120
3 59.42551 260 320 70
4 34.38484 200 158 110
5 33.17409 125 30 110
6 49.12025 200 125 90
7 53.31381 210 190 90
8 18.04285 220 35 120
9 50.76500 290 105 110
10 19.82357 210 45 120
11 40.40021 140 105 110
12 22.73645 180 55 110
13 41.44502 280 25 110
14 45.86332 290 35 100
15 35.78279 90 20 110
16 22.39651 180 65 110
17 64.53382 80 25 100
18 46.89564 220 30 110
19 44.33086 190 80 100
20 32.20758 125 30 110
21 31.43597 200 25 110
22 41.01549 240 190 120
23 28.02576 135 25 110
24 35.25244 45 40 100
25 23.80404 280 45 110
26 52.07690 140 85 100
27 53.37101 170 90 110
28 21.87129 220 45 120
29 31.07222 250 90 110
30 36.52368 170 60 110
31 39.24111 260 40 110
32 45.32807 150 95 100
33 26.73451 180 55 110
34 54.85092 0 95 100
35 40.10596 220 90 100
36 29.92429 190 40 120
37 30.45084 170 120 130
38 60.75611 0 15 50
39 63.00565 0 50 50
40 49.51187 135 110 100
41 50.82839 0 110 100
42 39.25920 210 240 120
43 3.70340 140 140 100
44 55.33314 0 110 90
45 41.99893 240 30 110
46 40.56016 290 35 110
47 68.23588 0 95 80
48 74.47295 0 140 90
49 72.80179 0 120 90
50 31.23005 70 40 110
51 53.13132 230 55 110
52 59.36399 15 90 90
53 38.83975 200 35 110
54 28.59278 190 230 140
55 46.65884 200 110 100
56 39.10617 250 60 110
57 27.75330 140 25 110
58 49.78744 230 115 100
59 51.59219 200 110 100
60 36.18756 200 60 110
> k <- length(x[1,])
> df <- as.data.frame(x)
> (mylm <- lm(df))
Call:
lm(formula = df)
Coefficients:
(Intercept) fat sugars calories
96.76424 -0.02523 0.04490 -0.52627
> (mysum <- summary(mylm))
Call:
lm(formula = df)
Residuals:
Min 1Q Median 3Q Max
-43.189 -6.823 0.255 6.338 21.292
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 96.76424 9.89233 9.782 1.01e-13 ***
fat -0.02523 0.01744 -1.447 0.1535
sugars 0.04490 0.02139 2.099 0.0403 *
calories -0.52627 0.09853 -5.341 1.74e-06 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 10.54 on 56 degrees of freedom
Multiple R-squared: 0.5044, Adjusted R-squared: 0.4779
F-statistic: 19 on 3 and 56 DF, p-value: 1.268e-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.046066736 0.092133472 0.95393326
[2,] 0.028337057 0.056674114 0.97166294
[3,] 0.345205389 0.690410778 0.65479461
[4,] 0.302621785 0.605243570 0.69737821
[5,] 0.217371514 0.434743028 0.78262849
[6,] 0.232782662 0.465565323 0.76721734
[7,] 0.225947384 0.451894768 0.77405262
[8,] 0.172338930 0.344677859 0.82766107
[9,] 0.115930730 0.231861460 0.88406927
[10,] 0.154809282 0.309618563 0.84519072
[11,] 0.375984256 0.751968513 0.62401574
[12,] 0.387752050 0.775504100 0.61224795
[13,] 0.305377721 0.610755442 0.69462228
[14,] 0.254132420 0.508264840 0.74586758
[15,] 0.202296485 0.404592970 0.79770352
[16,] 0.233695921 0.467391843 0.76630408
[17,] 0.209698482 0.419396965 0.79030152
[18,] 0.190382037 0.380764073 0.80961796
[19,] 0.194369146 0.388738292 0.80563085
[20,] 0.174065409 0.348130817 0.82593459
[21,] 0.242498175 0.484996350 0.75750182
[22,] 0.214037157 0.428074314 0.78596284
[23,] 0.171740788 0.343481575 0.82825921
[24,] 0.126518844 0.253037688 0.87348116
[25,] 0.095888629 0.191777258 0.90411137
[26,] 0.066397553 0.132795105 0.93360245
[27,] 0.064157632 0.128315265 0.93584237
[28,] 0.050819464 0.101638928 0.94918054
[29,] 0.034127078 0.068254155 0.96587292
[30,] 0.022296335 0.044592671 0.97770366
[31,] 0.015231902 0.030463803 0.98476810
[32,] 0.014639826 0.029279652 0.98536017
[33,] 0.015420780 0.030841560 0.98457922
[34,] 0.009735317 0.019470635 0.99026468
[35,] 0.005750382 0.011500763 0.99424962
[36,] 0.003137226 0.006274452 0.99686277
[37,] 0.944802400 0.110395201 0.05519760
[38,] 0.936498982 0.127002036 0.06350102
[39,] 0.919875740 0.160248520 0.08012426
[40,] 0.891715956 0.216568089 0.10828404
[41,] 0.862152848 0.275694304 0.13784715
[42,] 0.857123277 0.285753445 0.14287672
[43,] 0.912843870 0.174312259 0.08715613
[44,] 0.859527208 0.280945583 0.14047279
[45,] 0.986003347 0.027993305 0.01399665
[46,] 0.989806391 0.020387218 0.01019361
[47,] 0.989071278 0.021857444 0.01092872
> postscript(file="/var/wessaorg/rcomp/tmp/1nsvi1321882815.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/2mq651321882815.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/39m091321882815.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/4t3d81321882815.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/5gr3u1321882815.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.81584918 -5.31211385 -8.31000667 -6.53952649 -3.89479034 -0.84765677
7 8 9 10 11 12
0.67954276 -11.59141826 14.49077241 -10.51197522 0.34207479 -14.06755275
13 14 15 16 17 18
8.51067977 7.46957573 -1.71998606 -14.85650434 21.29161846 12.22323322
19 20 21 22 23 24
1.39392033 -4.86130234 -3.51644923 4.92595533 -8.56634967 -9.54619274
25 26 27 28 29 30
-10.02833342 7.65415016 14.74318280 -8.21199585 -5.53752420 -0.75708542
31 32 33 34 35 36
5.12872613 0.70856894 -10.06948375 6.44750632 -2.52320814 -0.69127468
37 38 39 40 41 42
1.00126466 -10.36840402 -9.69043609 3.84045050 1.75145344 0.16778626
43 44 45 46 47 48
-43.18894909 0.99355262 7.83104297 7.42906155 9.30717269 18.78630385
49 50 51 52 53 54
18.01317903 -7.67528099 17.58862712 6.30083137 3.43830518 -0.02882627
55 56 57 58 59 60
2.62711293 3.84348857 -8.71268348 6.28798576 7.56046193 -0.33642830
> postscript(file="/var/wessaorg/rcomp/tmp/6mi5i1321882815.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.81584918 NA
1 -5.31211385 -0.81584918
2 -8.31000667 -5.31211385
3 -6.53952649 -8.31000667
4 -3.89479034 -6.53952649
5 -0.84765677 -3.89479034
6 0.67954276 -0.84765677
7 -11.59141826 0.67954276
8 14.49077241 -11.59141826
9 -10.51197522 14.49077241
10 0.34207479 -10.51197522
11 -14.06755275 0.34207479
12 8.51067977 -14.06755275
13 7.46957573 8.51067977
14 -1.71998606 7.46957573
15 -14.85650434 -1.71998606
16 21.29161846 -14.85650434
17 12.22323322 21.29161846
18 1.39392033 12.22323322
19 -4.86130234 1.39392033
20 -3.51644923 -4.86130234
21 4.92595533 -3.51644923
22 -8.56634967 4.92595533
23 -9.54619274 -8.56634967
24 -10.02833342 -9.54619274
25 7.65415016 -10.02833342
26 14.74318280 7.65415016
27 -8.21199585 14.74318280
28 -5.53752420 -8.21199585
29 -0.75708542 -5.53752420
30 5.12872613 -0.75708542
31 0.70856894 5.12872613
32 -10.06948375 0.70856894
33 6.44750632 -10.06948375
34 -2.52320814 6.44750632
35 -0.69127468 -2.52320814
36 1.00126466 -0.69127468
37 -10.36840402 1.00126466
38 -9.69043609 -10.36840402
39 3.84045050 -9.69043609
40 1.75145344 3.84045050
41 0.16778626 1.75145344
42 -43.18894909 0.16778626
43 0.99355262 -43.18894909
44 7.83104297 0.99355262
45 7.42906155 7.83104297
46 9.30717269 7.42906155
47 18.78630385 9.30717269
48 18.01317903 18.78630385
49 -7.67528099 18.01317903
50 17.58862712 -7.67528099
51 6.30083137 17.58862712
52 3.43830518 6.30083137
53 -0.02882627 3.43830518
54 2.62711293 -0.02882627
55 3.84348857 2.62711293
56 -8.71268348 3.84348857
57 6.28798576 -8.71268348
58 7.56046193 6.28798576
59 -0.33642830 7.56046193
60 NA -0.33642830
> dum1 <- dum[2:length(myerror),]
> dum1
lag(myerror, k = 1) myerror
[1,] -5.31211385 -0.81584918
[2,] -8.31000667 -5.31211385
[3,] -6.53952649 -8.31000667
[4,] -3.89479034 -6.53952649
[5,] -0.84765677 -3.89479034
[6,] 0.67954276 -0.84765677
[7,] -11.59141826 0.67954276
[8,] 14.49077241 -11.59141826
[9,] -10.51197522 14.49077241
[10,] 0.34207479 -10.51197522
[11,] -14.06755275 0.34207479
[12,] 8.51067977 -14.06755275
[13,] 7.46957573 8.51067977
[14,] -1.71998606 7.46957573
[15,] -14.85650434 -1.71998606
[16,] 21.29161846 -14.85650434
[17,] 12.22323322 21.29161846
[18,] 1.39392033 12.22323322
[19,] -4.86130234 1.39392033
[20,] -3.51644923 -4.86130234
[21,] 4.92595533 -3.51644923
[22,] -8.56634967 4.92595533
[23,] -9.54619274 -8.56634967
[24,] -10.02833342 -9.54619274
[25,] 7.65415016 -10.02833342
[26,] 14.74318280 7.65415016
[27,] -8.21199585 14.74318280
[28,] -5.53752420 -8.21199585
[29,] -0.75708542 -5.53752420
[30,] 5.12872613 -0.75708542
[31,] 0.70856894 5.12872613
[32,] -10.06948375 0.70856894
[33,] 6.44750632 -10.06948375
[34,] -2.52320814 6.44750632
[35,] -0.69127468 -2.52320814
[36,] 1.00126466 -0.69127468
[37,] -10.36840402 1.00126466
[38,] -9.69043609 -10.36840402
[39,] 3.84045050 -9.69043609
[40,] 1.75145344 3.84045050
[41,] 0.16778626 1.75145344
[42,] -43.18894909 0.16778626
[43,] 0.99355262 -43.18894909
[44,] 7.83104297 0.99355262
[45,] 7.42906155 7.83104297
[46,] 9.30717269 7.42906155
[47,] 18.78630385 9.30717269
[48,] 18.01317903 18.78630385
[49,] -7.67528099 18.01317903
[50,] 17.58862712 -7.67528099
[51,] 6.30083137 17.58862712
[52,] 3.43830518 6.30083137
[53,] -0.02882627 3.43830518
[54,] 2.62711293 -0.02882627
[55,] 3.84348857 2.62711293
[56,] -8.71268348 3.84348857
[57,] 6.28798576 -8.71268348
[58,] 7.56046193 6.28798576
[59,] -0.33642830 7.56046193
> z <- as.data.frame(dum1)
> z
lag(myerror, k = 1) myerror
1 -5.31211385 -0.81584918
2 -8.31000667 -5.31211385
3 -6.53952649 -8.31000667
4 -3.89479034 -6.53952649
5 -0.84765677 -3.89479034
6 0.67954276 -0.84765677
7 -11.59141826 0.67954276
8 14.49077241 -11.59141826
9 -10.51197522 14.49077241
10 0.34207479 -10.51197522
11 -14.06755275 0.34207479
12 8.51067977 -14.06755275
13 7.46957573 8.51067977
14 -1.71998606 7.46957573
15 -14.85650434 -1.71998606
16 21.29161846 -14.85650434
17 12.22323322 21.29161846
18 1.39392033 12.22323322
19 -4.86130234 1.39392033
20 -3.51644923 -4.86130234
21 4.92595533 -3.51644923
22 -8.56634967 4.92595533
23 -9.54619274 -8.56634967
24 -10.02833342 -9.54619274
25 7.65415016 -10.02833342
26 14.74318280 7.65415016
27 -8.21199585 14.74318280
28 -5.53752420 -8.21199585
29 -0.75708542 -5.53752420
30 5.12872613 -0.75708542
31 0.70856894 5.12872613
32 -10.06948375 0.70856894
33 6.44750632 -10.06948375
34 -2.52320814 6.44750632
35 -0.69127468 -2.52320814
36 1.00126466 -0.69127468
37 -10.36840402 1.00126466
38 -9.69043609 -10.36840402
39 3.84045050 -9.69043609
40 1.75145344 3.84045050
41 0.16778626 1.75145344
42 -43.18894909 0.16778626
43 0.99355262 -43.18894909
44 7.83104297 0.99355262
45 7.42906155 7.83104297
46 9.30717269 7.42906155
47 18.78630385 9.30717269
48 18.01317903 18.78630385
49 -7.67528099 18.01317903
50 17.58862712 -7.67528099
51 6.30083137 17.58862712
52 3.43830518 6.30083137
53 -0.02882627 3.43830518
54 2.62711293 -0.02882627
55 3.84348857 2.62711293
56 -8.71268348 3.84348857
57 6.28798576 -8.71268348
58 7.56046193 6.28798576
59 -0.33642830 7.56046193
> 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/7cn1g1321882815.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/8cxw01321882815.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/9prjc1321882815.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/10iphs1321882815.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/11bi141321882815.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/12lfd21321882815.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/136gws1321882815.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/14z2bc1321882815.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/15qds21321882815.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/16ap3z1321882815.tab")
+ }
>
> try(system("convert tmp/1nsvi1321882815.ps tmp/1nsvi1321882815.png",intern=TRUE))
character(0)
> try(system("convert tmp/2mq651321882815.ps tmp/2mq651321882815.png",intern=TRUE))
character(0)
> try(system("convert tmp/39m091321882815.ps tmp/39m091321882815.png",intern=TRUE))
character(0)
> try(system("convert tmp/4t3d81321882815.ps tmp/4t3d81321882815.png",intern=TRUE))
character(0)
> try(system("convert tmp/5gr3u1321882815.ps tmp/5gr3u1321882815.png",intern=TRUE))
character(0)
> try(system("convert tmp/6mi5i1321882815.ps tmp/6mi5i1321882815.png",intern=TRUE))
character(0)
> try(system("convert tmp/7cn1g1321882815.ps tmp/7cn1g1321882815.png",intern=TRUE))
character(0)
> try(system("convert tmp/8cxw01321882815.ps tmp/8cxw01321882815.png",intern=TRUE))
character(0)
> try(system("convert tmp/9prjc1321882815.ps tmp/9prjc1321882815.png",intern=TRUE))
character(0)
> try(system("convert tmp/10iphs1321882815.ps tmp/10iphs1321882815.png",intern=TRUE))
character(0)
>
>
> proc.time()
user system elapsed
3.232 0.508 3.749