R version 2.8.0 (2008-10-20)
Copyright (C) 2008 The R Foundation for Statistical Computing
ISBN 3-900051-07-0
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.
Natural language support but running in an English locale
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(162556
+ ,807
+ ,213118
+ ,6282154
+ ,29790
+ ,444
+ ,81767
+ ,4321023
+ ,87550
+ ,412
+ ,153198
+ ,4111912
+ ,84738
+ ,428
+ ,-26007
+ ,223193
+ ,54660
+ ,315
+ ,126942
+ ,1491348
+ ,42634
+ ,168
+ ,157214
+ ,1629616
+ ,40949
+ ,263
+ ,129352
+ ,1398893
+ ,45187
+ ,267
+ ,234817
+ ,1926517
+ ,37704
+ ,228
+ ,60448
+ ,983660
+ ,16275
+ ,129
+ ,47818
+ ,1443586
+ ,25830
+ ,104
+ ,245546
+ ,1073089
+ ,12679
+ ,122
+ ,48020
+ ,984885
+ ,18014
+ ,393
+ ,-1710
+ ,1405225
+ ,43556
+ ,190
+ ,32648
+ ,227132
+ ,24811
+ ,280
+ ,95350
+ ,929118
+ ,6575
+ ,63
+ ,151352
+ ,1071292
+ ,7123
+ ,102
+ ,288170
+ ,638830
+ ,21950
+ ,265
+ ,114337
+ ,856956
+ ,37597
+ ,234
+ ,37884
+ ,992426
+ ,17821
+ ,277
+ ,122844
+ ,444477
+ ,12988
+ ,73
+ ,82340
+ ,857217
+ ,22330
+ ,67
+ ,79801
+ ,711969
+ ,13326
+ ,103
+ ,165548
+ ,702380
+ ,16189
+ ,290
+ ,116384
+ ,358589
+ ,7146
+ ,83
+ ,134028
+ ,297978
+ ,15824
+ ,56
+ ,63838
+ ,585715
+ ,27664
+ ,236
+ ,74996
+ ,657954
+ ,11920
+ ,73
+ ,31080
+ ,209458
+ ,8568
+ ,34
+ ,32168
+ ,786690
+ ,14416
+ ,139
+ ,49857
+ ,439798
+ ,3369
+ ,26
+ ,87161
+ ,688779
+ ,11819
+ ,70
+ ,106113
+ ,574339
+ ,6984
+ ,40
+ ,80570
+ ,741409
+ ,4519
+ ,42
+ ,102129
+ ,597793
+ ,2220
+ ,12
+ ,301670
+ ,644190
+ ,18562
+ ,211
+ ,102313
+ ,377934
+ ,10327
+ ,74
+ ,88577
+ ,640273
+ ,5336
+ ,80
+ ,112477
+ ,697458
+ ,2365
+ ,83
+ ,191778
+ ,550608
+ ,4069
+ ,131
+ ,79804
+ ,207393
+ ,8636
+ ,203
+ ,128294
+ ,301607
+ ,13718
+ ,56
+ ,96448
+ ,345783
+ ,4525
+ ,89
+ ,93811
+ ,501749
+ ,6869
+ ,88
+ ,117520
+ ,379983
+ ,4628
+ ,39
+ ,69159
+ ,387475
+ ,3689
+ ,25
+ ,101792
+ ,377305
+ ,4891
+ ,49
+ ,210568
+ ,370837
+ ,7489
+ ,149
+ ,136996
+ ,430866
+ ,4901
+ ,58
+ ,121920
+ ,469107
+ ,2284
+ ,41
+ ,76403
+ ,194493
+ ,3160
+ ,90
+ ,108094
+ ,530670
+ ,4150
+ ,136
+ ,134759
+ ,518365
+ ,7285
+ ,97
+ ,188873
+ ,491303
+ ,1134
+ ,63
+ ,146216
+ ,527021
+ ,4658
+ ,114
+ ,156608
+ ,233773
+ ,2384
+ ,77
+ ,61348
+ ,405972
+ ,3748
+ ,6
+ ,50350
+ ,652925
+ ,5371
+ ,47
+ ,87720
+ ,446211
+ ,1285
+ ,51
+ ,99489
+ ,341340
+ ,9327
+ ,85
+ ,87419
+ ,387699)
+ ,dim=c(4
+ ,60)
+ ,dimnames=list(c('Costs'
+ ,'Orders'
+ ,'Dividends'
+ ,'Wealth')
+ ,1:60))
> y <- array(NA,dim=c(4,60),dimnames=list(c('Costs','Orders','Dividends','Wealth'),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 = '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
Attaching package: 'zoo'
The following object(s) are masked from package:base :
as.Date.numeric
> n25 <- 25 #minimum number of obs. for Goldfeld-Quandt test
> par1 <- as.numeric(par1)
> x <- t(y)
> k <- length(x[1,])
> n <- length(x[,1])
> x1 <- cbind(x[,par1], x[,1:k!=par1])
> mycolnames <- c(colnames(x)[par1], colnames(x)[1:k!=par1])
> colnames(x1) <- mycolnames #colnames(x)[par1]
> x <- x1
> if (par3 == 'First Differences'){
+ x2 <- array(0, dim=c(n-1,k), dimnames=list(1:(n-1), paste('(1-B)',colnames(x),sep='')))
+ for (i in 1:n-1) {
+ for (j in 1:k) {
+ x2[i,j] <- x[i+1,j] - x[i,j]
+ }
+ }
+ x <- x2
+ }
> if (par2 == 'Include Monthly Dummies'){
+ x2 <- array(0, dim=c(n,11), dimnames=list(1:n, paste('M', seq(1:11), sep ='')))
+ for (i in 1:11){
+ x2[seq(i,n,12),i] <- 1
+ }
+ x <- cbind(x, x2)
+ }
> if (par2 == 'Include Quarterly Dummies'){
+ x2 <- array(0, dim=c(n,3), dimnames=list(1:n, paste('Q', seq(1:3), sep ='')))
+ for (i in 1:3){
+ x2[seq(i,n,4),i] <- 1
+ }
+ x <- cbind(x, x2)
+ }
> k <- length(x[1,])
> if (par3 == 'Linear Trend'){
+ x <- cbind(x, c(1:n))
+ colnames(x)[k+1] <- 't'
+ }
> x
Wealth Costs Orders Dividends t
1 6282154 162556 807 213118 1
2 4321023 29790 444 81767 2
3 4111912 87550 412 153198 3
4 223193 84738 428 -26007 4
5 1491348 54660 315 126942 5
6 1629616 42634 168 157214 6
7 1398893 40949 263 129352 7
8 1926517 45187 267 234817 8
9 983660 37704 228 60448 9
10 1443586 16275 129 47818 10
11 1073089 25830 104 245546 11
12 984885 12679 122 48020 12
13 1405225 18014 393 -1710 13
14 227132 43556 190 32648 14
15 929118 24811 280 95350 15
16 1071292 6575 63 151352 16
17 638830 7123 102 288170 17
18 856956 21950 265 114337 18
19 992426 37597 234 37884 19
20 444477 17821 277 122844 20
21 857217 12988 73 82340 21
22 711969 22330 67 79801 22
23 702380 13326 103 165548 23
24 358589 16189 290 116384 24
25 297978 7146 83 134028 25
26 585715 15824 56 63838 26
27 657954 27664 236 74996 27
28 209458 11920 73 31080 28
29 786690 8568 34 32168 29
30 439798 14416 139 49857 30
31 688779 3369 26 87161 31
32 574339 11819 70 106113 32
33 741409 6984 40 80570 33
34 597793 4519 42 102129 34
35 644190 2220 12 301670 35
36 377934 18562 211 102313 36
37 640273 10327 74 88577 37
38 697458 5336 80 112477 38
39 550608 2365 83 191778 39
40 207393 4069 131 79804 40
41 301607 8636 203 128294 41
42 345783 13718 56 96448 42
43 501749 4525 89 93811 43
44 379983 6869 88 117520 44
45 387475 4628 39 69159 45
46 377305 3689 25 101792 46
47 370837 4891 49 210568 47
48 430866 7489 149 136996 48
49 469107 4901 58 121920 49
50 194493 2284 41 76403 50
51 530670 3160 90 108094 51
52 518365 4150 136 134759 52
53 491303 7285 97 188873 53
54 527021 1134 63 146216 54
55 233773 4658 114 156608 55
56 405972 2384 77 61348 56
57 652925 3748 6 50350 57
58 446211 5371 47 87720 58
59 341340 1285 51 99489 59
60 387699 9327 85 87419 60
> k <- length(x[1,])
> df <- as.data.frame(x)
> (mylm <- lm(df))
Call:
lm(formula = df)
Coefficients:
(Intercept) Costs Orders Dividends t
-42000.819 15.921 2669.759 2.988 -4099.567
> (mysum <- summary(mylm))
Call:
lm(formula = df)
Residuals:
Min 1Q Median 3Q Max
-2132506 -264617 -9646 234791 2467210
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) -42000.819 314242.402 -0.134 0.8942
Costs 15.921 6.324 2.518 0.0148 *
Orders 2669.759 1168.612 2.285 0.0262 *
Dividends 2.988 1.272 2.350 0.0224 *
t -4099.567 6197.912 -0.661 0.5111
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Residual standard error: 618100 on 55 degrees of freedom
Multiple R-squared: 0.6652, Adjusted R-squared: 0.6408
F-statistic: 27.32 on 4 and 55 DF, p-value: 1.649e-12
> 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.9971405 5.719052e-03 2.859526e-03
[2,] 1.0000000 3.113197e-08 1.556598e-08
[3,] 1.0000000 1.621720e-12 8.108600e-13
[4,] 1.0000000 2.109307e-12 1.054653e-12
[5,] 1.0000000 1.542027e-12 7.710134e-13
[6,] 1.0000000 6.536047e-14 3.268024e-14
[7,] 1.0000000 3.739284e-15 1.869642e-15
[8,] 1.0000000 5.281569e-15 2.640784e-15
[9,] 1.0000000 5.954250e-16 2.977125e-16
[10,] 1.0000000 9.262105e-16 4.631052e-16
[11,] 1.0000000 1.122745e-15 5.613725e-16
[12,] 1.0000000 2.622950e-16 1.311475e-16
[13,] 1.0000000 8.998878e-16 4.499439e-16
[14,] 1.0000000 1.453233e-16 7.266165e-17
[15,] 1.0000000 2.540725e-16 1.270362e-16
[16,] 1.0000000 1.117470e-15 5.587351e-16
[17,] 1.0000000 4.387580e-15 2.193790e-15
[18,] 1.0000000 5.192560e-15 2.596280e-15
[19,] 1.0000000 1.237282e-14 6.186412e-15
[20,] 1.0000000 2.018906e-14 1.009453e-14
[21,] 1.0000000 3.372440e-15 1.686220e-15
[22,] 1.0000000 3.662202e-15 1.831101e-15
[23,] 1.0000000 2.211629e-14 1.105814e-14
[24,] 1.0000000 8.027079e-14 4.013540e-14
[25,] 1.0000000 4.789234e-13 2.394617e-13
[26,] 1.0000000 9.054738e-13 4.527369e-13
[27,] 1.0000000 4.471446e-12 2.235723e-12
[28,] 1.0000000 2.623230e-11 1.311615e-11
[29,] 1.0000000 1.402929e-10 7.014645e-11
[30,] 1.0000000 2.850716e-10 1.425358e-10
[31,] 1.0000000 1.459395e-10 7.296977e-11
[32,] 1.0000000 4.574262e-10 2.287131e-10
[33,] 1.0000000 1.522656e-09 7.613279e-10
[34,] 1.0000000 7.781210e-09 3.890605e-09
[35,] 1.0000000 4.384447e-08 2.192224e-08
[36,] 0.9999999 1.910510e-07 9.552551e-08
[37,] 0.9999994 1.105004e-06 5.525022e-07
[38,] 0.9999971 5.767229e-06 2.883615e-06
[39,] 0.9999861 2.786107e-05 1.393054e-05
[40,] 0.9999398 1.204167e-04 6.020836e-05
[41,] 0.9997154 5.692546e-04 2.846273e-04
[42,] 0.9987469 2.506271e-03 1.253136e-03
[43,] 0.9999184 1.632386e-04 8.161931e-05
[44,] 0.9995603 8.794783e-04 4.397392e-04
[45,] 0.9987813 2.437348e-03 1.218674e-03
> postscript(file="/var/www/html/freestat/rcomp/tmp/11d061291303119.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/html/freestat/rcomp/tmp/21d061291303119.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/html/freestat/rcomp/tmp/3umz91291303119.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/html/freestat/rcomp/tmp/4umz91291303119.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/html/freestat/rcomp/tmp/5umz91291303119.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
948778.42 2467210.05 1214552.93 -2132505.81 -536728.28 99106.11
7 8 9 10 11 12
-271058.65 -132643.00 -327084.76 780170.07 -262469.86 405006.60
13 14 15 16 17 18
169604.95 -971763.93 -394876.10 453732.80 -496318.73 -425879.11
19 20 21 22 23 24
-224211.64 -821877.07 337577.13 71295.72 -143179.52 -880785.82
25 26 27 28 29 30
-293403.37 142093.71 -483977.52 -111304.12 624265.72 -144819.10
31 32 33 34 35 36
474356.12 55376.71 459947.83 289915.14 -139166.20 -597069.41
37 38 39 40 41 42
207285.98 260597.52 -79828.97 -239619.10 -551140.88 -96160.68
43 44 45 46 47 48
130049.14 -93115.45 229488.13 178229.73 -232398.15 -256759.54
49 50 51 52 53 54
114784.46 67337.71 168148.89 -58308.86 -188769.62 167221.78
55 56 57 58 59 60
-345245.40 250698.43 702451.50 252866.24 171302.55 39016.55
> postscript(file="/var/www/html/freestat/rcomp/tmp/65dgc1291303119.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 948778.42 NA
1 2467210.05 948778.42
2 1214552.93 2467210.05
3 -2132505.81 1214552.93
4 -536728.28 -2132505.81
5 99106.11 -536728.28
6 -271058.65 99106.11
7 -132643.00 -271058.65
8 -327084.76 -132643.00
9 780170.07 -327084.76
10 -262469.86 780170.07
11 405006.60 -262469.86
12 169604.95 405006.60
13 -971763.93 169604.95
14 -394876.10 -971763.93
15 453732.80 -394876.10
16 -496318.73 453732.80
17 -425879.11 -496318.73
18 -224211.64 -425879.11
19 -821877.07 -224211.64
20 337577.13 -821877.07
21 71295.72 337577.13
22 -143179.52 71295.72
23 -880785.82 -143179.52
24 -293403.37 -880785.82
25 142093.71 -293403.37
26 -483977.52 142093.71
27 -111304.12 -483977.52
28 624265.72 -111304.12
29 -144819.10 624265.72
30 474356.12 -144819.10
31 55376.71 474356.12
32 459947.83 55376.71
33 289915.14 459947.83
34 -139166.20 289915.14
35 -597069.41 -139166.20
36 207285.98 -597069.41
37 260597.52 207285.98
38 -79828.97 260597.52
39 -239619.10 -79828.97
40 -551140.88 -239619.10
41 -96160.68 -551140.88
42 130049.14 -96160.68
43 -93115.45 130049.14
44 229488.13 -93115.45
45 178229.73 229488.13
46 -232398.15 178229.73
47 -256759.54 -232398.15
48 114784.46 -256759.54
49 67337.71 114784.46
50 168148.89 67337.71
51 -58308.86 168148.89
52 -188769.62 -58308.86
53 167221.78 -188769.62
54 -345245.40 167221.78
55 250698.43 -345245.40
56 702451.50 250698.43
57 252866.24 702451.50
58 171302.55 252866.24
59 39016.55 171302.55
60 NA 39016.55
> dum1 <- dum[2:length(myerror),]
> dum1
lag(myerror, k = 1) myerror
[1,] 2467210.05 948778.42
[2,] 1214552.93 2467210.05
[3,] -2132505.81 1214552.93
[4,] -536728.28 -2132505.81
[5,] 99106.11 -536728.28
[6,] -271058.65 99106.11
[7,] -132643.00 -271058.65
[8,] -327084.76 -132643.00
[9,] 780170.07 -327084.76
[10,] -262469.86 780170.07
[11,] 405006.60 -262469.86
[12,] 169604.95 405006.60
[13,] -971763.93 169604.95
[14,] -394876.10 -971763.93
[15,] 453732.80 -394876.10
[16,] -496318.73 453732.80
[17,] -425879.11 -496318.73
[18,] -224211.64 -425879.11
[19,] -821877.07 -224211.64
[20,] 337577.13 -821877.07
[21,] 71295.72 337577.13
[22,] -143179.52 71295.72
[23,] -880785.82 -143179.52
[24,] -293403.37 -880785.82
[25,] 142093.71 -293403.37
[26,] -483977.52 142093.71
[27,] -111304.12 -483977.52
[28,] 624265.72 -111304.12
[29,] -144819.10 624265.72
[30,] 474356.12 -144819.10
[31,] 55376.71 474356.12
[32,] 459947.83 55376.71
[33,] 289915.14 459947.83
[34,] -139166.20 289915.14
[35,] -597069.41 -139166.20
[36,] 207285.98 -597069.41
[37,] 260597.52 207285.98
[38,] -79828.97 260597.52
[39,] -239619.10 -79828.97
[40,] -551140.88 -239619.10
[41,] -96160.68 -551140.88
[42,] 130049.14 -96160.68
[43,] -93115.45 130049.14
[44,] 229488.13 -93115.45
[45,] 178229.73 229488.13
[46,] -232398.15 178229.73
[47,] -256759.54 -232398.15
[48,] 114784.46 -256759.54
[49,] 67337.71 114784.46
[50,] 168148.89 67337.71
[51,] -58308.86 168148.89
[52,] -188769.62 -58308.86
[53,] 167221.78 -188769.62
[54,] -345245.40 167221.78
[55,] 250698.43 -345245.40
[56,] 702451.50 250698.43
[57,] 252866.24 702451.50
[58,] 171302.55 252866.24
[59,] 39016.55 171302.55
> z <- as.data.frame(dum1)
> z
lag(myerror, k = 1) myerror
1 2467210.05 948778.42
2 1214552.93 2467210.05
3 -2132505.81 1214552.93
4 -536728.28 -2132505.81
5 99106.11 -536728.28
6 -271058.65 99106.11
7 -132643.00 -271058.65
8 -327084.76 -132643.00
9 780170.07 -327084.76
10 -262469.86 780170.07
11 405006.60 -262469.86
12 169604.95 405006.60
13 -971763.93 169604.95
14 -394876.10 -971763.93
15 453732.80 -394876.10
16 -496318.73 453732.80
17 -425879.11 -496318.73
18 -224211.64 -425879.11
19 -821877.07 -224211.64
20 337577.13 -821877.07
21 71295.72 337577.13
22 -143179.52 71295.72
23 -880785.82 -143179.52
24 -293403.37 -880785.82
25 142093.71 -293403.37
26 -483977.52 142093.71
27 -111304.12 -483977.52
28 624265.72 -111304.12
29 -144819.10 624265.72
30 474356.12 -144819.10
31 55376.71 474356.12
32 459947.83 55376.71
33 289915.14 459947.83
34 -139166.20 289915.14
35 -597069.41 -139166.20
36 207285.98 -597069.41
37 260597.52 207285.98
38 -79828.97 260597.52
39 -239619.10 -79828.97
40 -551140.88 -239619.10
41 -96160.68 -551140.88
42 130049.14 -96160.68
43 -93115.45 130049.14
44 229488.13 -93115.45
45 178229.73 229488.13
46 -232398.15 178229.73
47 -256759.54 -232398.15
48 114784.46 -256759.54
49 67337.71 114784.46
50 168148.89 67337.71
51 -58308.86 168148.89
52 -188769.62 -58308.86
53 167221.78 -188769.62
54 -345245.40 167221.78
55 250698.43 -345245.40
56 702451.50 250698.43
57 252866.24 702451.50
58 171302.55 252866.24
59 39016.55 171302.55
> 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/html/freestat/rcomp/tmp/75dgc1291303119.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/html/freestat/rcomp/tmp/8x5yx1291303119.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/html/freestat/rcomp/tmp/9x5yx1291303119.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/html/freestat/rcomp/tmp/10qwx01291303119.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/html/freestat/rcomp/createtable file can be downloaded at http://www.wessa.net/cretab
> load(file="/var/www/html/freestat/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/html/freestat/rcomp/tmp/11twd61291303119.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/html/freestat/rcomp/tmp/12xxcu1291303119.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/html/freestat/rcomp/tmp/13t7al1291303119.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/html/freestat/rcomp/tmp/14w7891291303119.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/html/freestat/rcomp/tmp/15087e1291303119.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/html/freestat/rcomp/tmp/1638521291303119.tab")
+ }
>
> try(system("convert tmp/11d061291303119.ps tmp/11d061291303119.png",intern=TRUE))
character(0)
> try(system("convert tmp/21d061291303119.ps tmp/21d061291303119.png",intern=TRUE))
character(0)
> try(system("convert tmp/3umz91291303119.ps tmp/3umz91291303119.png",intern=TRUE))
character(0)
> try(system("convert tmp/4umz91291303119.ps tmp/4umz91291303119.png",intern=TRUE))
character(0)
> try(system("convert tmp/5umz91291303119.ps tmp/5umz91291303119.png",intern=TRUE))
character(0)
> try(system("convert tmp/65dgc1291303119.ps tmp/65dgc1291303119.png",intern=TRUE))
character(0)
> try(system("convert tmp/75dgc1291303119.ps tmp/75dgc1291303119.png",intern=TRUE))
character(0)
> try(system("convert tmp/8x5yx1291303119.ps tmp/8x5yx1291303119.png",intern=TRUE))
character(0)
> try(system("convert tmp/9x5yx1291303119.ps tmp/9x5yx1291303119.png",intern=TRUE))
character(0)
> try(system("convert tmp/10qwx01291303119.ps tmp/10qwx01291303119.png",intern=TRUE))
character(0)
>
>
> proc.time()
user system elapsed
3.871 2.470 4.222