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 = '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
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
1 6282154 162556 807 213118
2 4321023 29790 444 81767
3 4111912 87550 412 153198
4 223193 84738 428 -26007
5 1491348 54660 315 126942
6 1629616 42634 168 157214
7 1398893 40949 263 129352
8 1926517 45187 267 234817
9 983660 37704 228 60448
10 1443586 16275 129 47818
11 1073089 25830 104 245546
12 984885 12679 122 48020
13 1405225 18014 393 -1710
14 227132 43556 190 32648
15 929118 24811 280 95350
16 1071292 6575 63 151352
17 638830 7123 102 288170
18 856956 21950 265 114337
19 992426 37597 234 37884
20 444477 17821 277 122844
21 857217 12988 73 82340
22 711969 22330 67 79801
23 702380 13326 103 165548
24 358589 16189 290 116384
25 297978 7146 83 134028
26 585715 15824 56 63838
27 657954 27664 236 74996
28 209458 11920 73 31080
29 786690 8568 34 32168
30 439798 14416 139 49857
31 688779 3369 26 87161
32 574339 11819 70 106113
33 741409 6984 40 80570
34 597793 4519 42 102129
35 644190 2220 12 301670
36 377934 18562 211 102313
37 640273 10327 74 88577
38 697458 5336 80 112477
39 550608 2365 83 191778
40 207393 4069 131 79804
41 301607 8636 203 128294
42 345783 13718 56 96448
43 501749 4525 89 93811
44 379983 6869 88 117520
45 387475 4628 39 69159
46 377305 3689 25 101792
47 370837 4891 49 210568
48 430866 7489 149 136996
49 469107 4901 58 121920
50 194493 2284 41 76403
51 530670 3160 90 108094
52 518365 4150 136 134759
53 491303 7285 97 188873
54 527021 1134 63 146216
55 233773 4658 114 156608
56 405972 2384 77 61348
57 652925 3748 6 50350
58 446211 5371 47 87720
59 341340 1285 51 99489
60 387699 9327 85 87419
> k <- length(x[1,])
> df <- as.data.frame(x)
> (mylm <- lm(df))
Call:
lm(formula = df)
Coefficients:
(Intercept) Costs Orders Dividends
-2.081e+05 1.699e+01 2.823e+03 2.968e+00
> (mysum <- summary(mylm))
Call:
lm(formula = df)
Residuals:
Min 1Q Median 3Q Max
-2139838 -262080 -72691 201243 2526719
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) -2.081e+05 1.879e+05 -1.108 0.27268
Costs 1.699e+01 6.082e+00 2.794 0.00711 **
Orders 2.823e+03 1.140e+03 2.478 0.01627 *
Dividends 2.968e+00 1.265e+00 2.346 0.02253 *
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Residual standard error: 615000 on 56 degrees of freedom
Multiple R-squared: 0.6625, Adjusted R-squared: 0.6445
F-statistic: 36.65 on 3 and 56 DF, p-value: 3.067e-13
> 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.9999863 2.742075e-05 1.371038e-05
[2,] 1.0000000 5.002820e-08 2.501410e-08
[3,] 0.9999999 1.439989e-07 7.199946e-08
[4,] 1.0000000 3.427399e-10 1.713700e-10
[5,] 1.0000000 2.929078e-10 1.464539e-10
[6,] 1.0000000 1.978984e-10 9.894921e-11
[7,] 1.0000000 3.210496e-13 1.605248e-13
[8,] 1.0000000 1.234901e-14 6.174505e-15
[9,] 1.0000000 2.823086e-15 1.411543e-15
[10,] 1.0000000 1.749288e-16 8.746441e-17
[11,] 1.0000000 5.947393e-17 2.973697e-17
[12,] 1.0000000 1.846476e-17 9.232382e-18
[13,] 1.0000000 4.172260e-17 2.086130e-17
[14,] 1.0000000 2.852901e-17 1.426450e-17
[15,] 1.0000000 1.628021e-17 8.140105e-18
[16,] 1.0000000 7.383644e-17 3.691822e-17
[17,] 1.0000000 2.810869e-16 1.405435e-16
[18,] 1.0000000 3.983363e-16 1.991681e-16
[19,] 1.0000000 1.202857e-15 6.014287e-16
[20,] 1.0000000 5.714381e-15 2.857190e-15
[21,] 1.0000000 1.069248e-14 5.346239e-15
[22,] 1.0000000 9.423299e-15 4.711650e-15
[23,] 1.0000000 7.006857e-15 3.503429e-15
[24,] 1.0000000 4.350563e-14 2.175282e-14
[25,] 1.0000000 8.971372e-14 4.485686e-14
[26,] 1.0000000 4.848342e-13 2.424171e-13
[27,] 1.0000000 4.426645e-13 2.213322e-13
[28,] 1.0000000 1.610125e-12 8.050627e-13
[29,] 1.0000000 8.273186e-12 4.136593e-12
[30,] 1.0000000 3.820888e-11 1.910444e-11
[31,] 1.0000000 6.373628e-11 3.186814e-11
[32,] 1.0000000 3.538023e-11 1.769012e-11
[33,] 1.0000000 1.573507e-10 7.867536e-11
[34,] 1.0000000 3.260447e-10 1.630224e-10
[35,] 1.0000000 1.387075e-09 6.935374e-10
[36,] 1.0000000 7.805514e-09 3.902757e-09
[37,] 1.0000000 4.012213e-08 2.006106e-08
[38,] 0.9999999 2.256361e-07 1.128181e-07
[39,] 0.9999994 1.191578e-06 5.957891e-07
[40,] 0.9999970 5.911162e-06 2.955581e-06
[41,] 0.9999868 2.639251e-05 1.319626e-05
[42,] 0.9999364 1.272981e-04 6.364903e-05
[43,] 0.9997002 5.996028e-04 2.998014e-04
[44,] 0.9997112 5.775610e-04 2.887805e-04
[45,] 0.9988228 2.354337e-03 1.177168e-03
[46,] 0.9989856 2.028834e-03 1.014417e-03
[47,] 0.9949431 1.011372e-02 5.056860e-03
> postscript(file="/var/www/html/freestat/rcomp/tmp/187ay1291121342.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/2jy911291121342.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/3jy911291121342.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/4jy911291121342.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/5cqqm1291121342.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
817033.66 2526719.11 1214426.67 -2139837.63 -495444.32 172369.65
7 8 9 10 11 12
-215244.52 -83929.10 -272026.30 869039.63 -180060.75 490610.26
13 14 15 16 17 18
202764.81 -938209.46 -357866.48 540651.88 -417277.78 -395410.43
19 20 21 22 23 24
-211416.67 -796850.50 394178.53 114653.80 -98044.79 -872536.36
25 26 27 28 29 30
-247417.36 177389.51 -492880.97 -83302.83 657770.54 -137442.07
31 32 33 34 35 36
507586.78 69082.45 478820.37 307463.83 -114567.66 -628716.19
37 38 39 40 41 42
201120.41 255249.25 -84931.91 -260306.95 -590887.02 -123534.05
43 44 45 46 47 48
103308.86 -125829.47 201611.38 150076.52 -267400.62 -315504.20
49 50 51 52 53 54
68378.12 21316.51 110213.57 -127922.63 -258748.88 104084.17
55 56 57 58 59 60
-423877.83 174137.44 631005.49 170049.85 88392.54 -62079.84
> postscript(file="/var/www/html/freestat/rcomp/tmp/6cqqm1291121342.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 817033.66 NA
1 2526719.11 817033.66
2 1214426.67 2526719.11
3 -2139837.63 1214426.67
4 -495444.32 -2139837.63
5 172369.65 -495444.32
6 -215244.52 172369.65
7 -83929.10 -215244.52
8 -272026.30 -83929.10
9 869039.63 -272026.30
10 -180060.75 869039.63
11 490610.26 -180060.75
12 202764.81 490610.26
13 -938209.46 202764.81
14 -357866.48 -938209.46
15 540651.88 -357866.48
16 -417277.78 540651.88
17 -395410.43 -417277.78
18 -211416.67 -395410.43
19 -796850.50 -211416.67
20 394178.53 -796850.50
21 114653.80 394178.53
22 -98044.79 114653.80
23 -872536.36 -98044.79
24 -247417.36 -872536.36
25 177389.51 -247417.36
26 -492880.97 177389.51
27 -83302.83 -492880.97
28 657770.54 -83302.83
29 -137442.07 657770.54
30 507586.78 -137442.07
31 69082.45 507586.78
32 478820.37 69082.45
33 307463.83 478820.37
34 -114567.66 307463.83
35 -628716.19 -114567.66
36 201120.41 -628716.19
37 255249.25 201120.41
38 -84931.91 255249.25
39 -260306.95 -84931.91
40 -590887.02 -260306.95
41 -123534.05 -590887.02
42 103308.86 -123534.05
43 -125829.47 103308.86
44 201611.38 -125829.47
45 150076.52 201611.38
46 -267400.62 150076.52
47 -315504.20 -267400.62
48 68378.12 -315504.20
49 21316.51 68378.12
50 110213.57 21316.51
51 -127922.63 110213.57
52 -258748.88 -127922.63
53 104084.17 -258748.88
54 -423877.83 104084.17
55 174137.44 -423877.83
56 631005.49 174137.44
57 170049.85 631005.49
58 88392.54 170049.85
59 -62079.84 88392.54
60 NA -62079.84
> dum1 <- dum[2:length(myerror),]
> dum1
lag(myerror, k = 1) myerror
[1,] 2526719.11 817033.66
[2,] 1214426.67 2526719.11
[3,] -2139837.63 1214426.67
[4,] -495444.32 -2139837.63
[5,] 172369.65 -495444.32
[6,] -215244.52 172369.65
[7,] -83929.10 -215244.52
[8,] -272026.30 -83929.10
[9,] 869039.63 -272026.30
[10,] -180060.75 869039.63
[11,] 490610.26 -180060.75
[12,] 202764.81 490610.26
[13,] -938209.46 202764.81
[14,] -357866.48 -938209.46
[15,] 540651.88 -357866.48
[16,] -417277.78 540651.88
[17,] -395410.43 -417277.78
[18,] -211416.67 -395410.43
[19,] -796850.50 -211416.67
[20,] 394178.53 -796850.50
[21,] 114653.80 394178.53
[22,] -98044.79 114653.80
[23,] -872536.36 -98044.79
[24,] -247417.36 -872536.36
[25,] 177389.51 -247417.36
[26,] -492880.97 177389.51
[27,] -83302.83 -492880.97
[28,] 657770.54 -83302.83
[29,] -137442.07 657770.54
[30,] 507586.78 -137442.07
[31,] 69082.45 507586.78
[32,] 478820.37 69082.45
[33,] 307463.83 478820.37
[34,] -114567.66 307463.83
[35,] -628716.19 -114567.66
[36,] 201120.41 -628716.19
[37,] 255249.25 201120.41
[38,] -84931.91 255249.25
[39,] -260306.95 -84931.91
[40,] -590887.02 -260306.95
[41,] -123534.05 -590887.02
[42,] 103308.86 -123534.05
[43,] -125829.47 103308.86
[44,] 201611.38 -125829.47
[45,] 150076.52 201611.38
[46,] -267400.62 150076.52
[47,] -315504.20 -267400.62
[48,] 68378.12 -315504.20
[49,] 21316.51 68378.12
[50,] 110213.57 21316.51
[51,] -127922.63 110213.57
[52,] -258748.88 -127922.63
[53,] 104084.17 -258748.88
[54,] -423877.83 104084.17
[55,] 174137.44 -423877.83
[56,] 631005.49 174137.44
[57,] 170049.85 631005.49
[58,] 88392.54 170049.85
[59,] -62079.84 88392.54
> z <- as.data.frame(dum1)
> z
lag(myerror, k = 1) myerror
1 2526719.11 817033.66
2 1214426.67 2526719.11
3 -2139837.63 1214426.67
4 -495444.32 -2139837.63
5 172369.65 -495444.32
6 -215244.52 172369.65
7 -83929.10 -215244.52
8 -272026.30 -83929.10
9 869039.63 -272026.30
10 -180060.75 869039.63
11 490610.26 -180060.75
12 202764.81 490610.26
13 -938209.46 202764.81
14 -357866.48 -938209.46
15 540651.88 -357866.48
16 -417277.78 540651.88
17 -395410.43 -417277.78
18 -211416.67 -395410.43
19 -796850.50 -211416.67
20 394178.53 -796850.50
21 114653.80 394178.53
22 -98044.79 114653.80
23 -872536.36 -98044.79
24 -247417.36 -872536.36
25 177389.51 -247417.36
26 -492880.97 177389.51
27 -83302.83 -492880.97
28 657770.54 -83302.83
29 -137442.07 657770.54
30 507586.78 -137442.07
31 69082.45 507586.78
32 478820.37 69082.45
33 307463.83 478820.37
34 -114567.66 307463.83
35 -628716.19 -114567.66
36 201120.41 -628716.19
37 255249.25 201120.41
38 -84931.91 255249.25
39 -260306.95 -84931.91
40 -590887.02 -260306.95
41 -123534.05 -590887.02
42 103308.86 -123534.05
43 -125829.47 103308.86
44 201611.38 -125829.47
45 150076.52 201611.38
46 -267400.62 150076.52
47 -315504.20 -267400.62
48 68378.12 -315504.20
49 21316.51 68378.12
50 110213.57 21316.51
51 -127922.63 110213.57
52 -258748.88 -127922.63
53 104084.17 -258748.88
54 -423877.83 104084.17
55 174137.44 -423877.83
56 631005.49 174137.44
57 170049.85 631005.49
58 88392.54 170049.85
59 -62079.84 88392.54
> 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/74hpp1291121342.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/84hpp1291121342.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/9x8pa1291121342.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/10x8pa1291121342.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/111r5x1291121342.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/12mr431291121342.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/13bsjf1291121342.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/1431i01291121342.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/15pkg61291121342.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/16lcef1291121342.tab")
+ }
>
> try(system("convert tmp/187ay1291121342.ps tmp/187ay1291121342.png",intern=TRUE))
character(0)
> try(system("convert tmp/2jy911291121342.ps tmp/2jy911291121342.png",intern=TRUE))
character(0)
> try(system("convert tmp/3jy911291121342.ps tmp/3jy911291121342.png",intern=TRUE))
character(0)
> try(system("convert tmp/4jy911291121342.ps tmp/4jy911291121342.png",intern=TRUE))
character(0)
> try(system("convert tmp/5cqqm1291121342.ps tmp/5cqqm1291121342.png",intern=TRUE))
character(0)
> try(system("convert tmp/6cqqm1291121342.ps tmp/6cqqm1291121342.png",intern=TRUE))
character(0)
> try(system("convert tmp/74hpp1291121342.ps tmp/74hpp1291121342.png",intern=TRUE))
character(0)
> try(system("convert tmp/84hpp1291121342.ps tmp/84hpp1291121342.png",intern=TRUE))
character(0)
> try(system("convert tmp/9x8pa1291121342.ps tmp/9x8pa1291121342.png",intern=TRUE))
character(0)
> try(system("convert tmp/10x8pa1291121342.ps tmp/10x8pa1291121342.png",intern=TRUE))
character(0)
>
>
> proc.time()
user system elapsed
3.868 2.479 4.190