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(2
+ ,210907
+ ,79
+ ,94
+ ,4
+ ,179321
+ ,108
+ ,103
+ ,-4
+ ,173326
+ ,86
+ ,148
+ ,4
+ ,133131
+ ,44
+ ,90
+ ,4
+ ,258873
+ ,104
+ ,124
+ ,-1
+ ,230964
+ ,102
+ ,115
+ ,1
+ ,344297
+ ,80
+ ,108
+ ,3
+ ,174415
+ ,73
+ ,114
+ ,-1
+ ,223632
+ ,105
+ ,120
+ ,4
+ ,294424
+ ,107
+ ,124
+ ,3
+ ,325107
+ ,84
+ ,126
+ ,1
+ ,106408
+ ,33
+ ,37
+ ,-2
+ ,265769
+ ,96
+ ,120
+ ,-3
+ ,269651
+ ,106
+ ,93
+ ,-4
+ ,149112
+ ,56
+ ,95
+ ,2
+ ,152871
+ ,59
+ ,90
+ ,2
+ ,362301
+ ,76
+ ,110
+ ,-4
+ ,183167
+ ,91
+ ,138
+ ,3
+ ,277965
+ ,115
+ ,133
+ ,2
+ ,218946
+ ,76
+ ,96
+ ,2
+ ,244052
+ ,101
+ ,164
+ ,5
+ ,233328
+ ,92
+ ,102
+ ,-2
+ ,206161
+ ,75
+ ,99
+ ,-2
+ ,207176
+ ,56
+ ,114
+ ,-3
+ ,196553
+ ,41
+ ,99
+ ,2
+ ,143246
+ ,67
+ ,104
+ ,2
+ ,182192
+ ,77
+ ,138
+ ,2
+ ,194979
+ ,66
+ ,151
+ ,4
+ ,143756
+ ,105
+ ,120
+ ,4
+ ,275541
+ ,116
+ ,115
+ ,2
+ ,152299
+ ,62
+ ,98
+ ,2
+ ,193339
+ ,100
+ ,71
+ ,-4
+ ,130585
+ ,67
+ ,107
+ ,3
+ ,112611
+ ,46
+ ,73
+ ,3
+ ,148446
+ ,135
+ ,129
+ ,2
+ ,182079
+ ,124
+ ,118
+ ,-1
+ ,243060
+ ,58
+ ,104
+ ,-3
+ ,162765
+ ,68
+ ,107
+ ,1
+ ,225060
+ ,93
+ ,139
+ ,-3
+ ,133328
+ ,56
+ ,56
+ ,3
+ ,100750
+ ,83
+ ,93
+ ,3
+ ,132487
+ ,71
+ ,98
+ ,-3
+ ,317394
+ ,116
+ ,82
+ ,-4
+ ,184510
+ ,64
+ ,140
+ ,2
+ ,128423
+ ,32
+ ,120
+ ,-1
+ ,97839
+ ,25
+ ,66
+ ,3
+ ,172494
+ ,46
+ ,139
+ ,2
+ ,229242
+ ,63
+ ,119
+ ,5
+ ,351619
+ ,95
+ ,141
+ ,2
+ ,324598
+ ,113
+ ,133
+ ,-2
+ ,195838
+ ,111
+ ,98
+ ,3
+ ,199476
+ ,87
+ ,105
+ ,-2
+ ,92499
+ ,25
+ ,55
+ ,6
+ ,181633
+ ,47
+ ,73
+ ,-3
+ ,271856
+ ,109
+ ,86
+ ,3
+ ,95227
+ ,37
+ ,48
+ ,-2
+ ,118612
+ ,54
+ ,43
+ ,1
+ ,65475
+ ,16
+ ,46
+ ,2
+ ,121848
+ ,37
+ ,52
+ ,2
+ ,76302
+ ,29
+ ,68
+ ,-3
+ ,98104
+ ,55
+ ,47
+ ,-2
+ ,30989
+ ,5
+ ,41
+ ,1
+ ,31774
+ ,0
+ ,47
+ ,-4
+ ,150580
+ ,27
+ ,71
+ ,1
+ ,59382
+ ,29
+ ,24)
+ ,dim=c(4
+ ,65)
+ ,dimnames=list(c('Tot'
+ ,'Time'
+ ,'Comp'
+ ,'LFB')
+ ,1:65))
> y <- array(NA,dim=c(4,65),dimnames=list(c('Tot','Time','Comp','LFB'),1:65))
> for (i in 1:dim(x)[1])
+ {
+ for (j in 1:dim(x)[2])
+ {
+ y[i,j] <- as.numeric(x[i,j])
+ }
+ }
> par3 = 'No Linear Trend'
> par2 = 'Do not include Seasonal Dummies'
> par1 = '1'
> 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
Tot Time Comp LFB
1 2 210907 79 94
2 4 179321 108 103
3 -4 173326 86 148
4 4 133131 44 90
5 4 258873 104 124
6 -1 230964 102 115
7 1 344297 80 108
8 3 174415 73 114
9 -1 223632 105 120
10 4 294424 107 124
11 3 325107 84 126
12 1 106408 33 37
13 -2 265769 96 120
14 -3 269651 106 93
15 -4 149112 56 95
16 2 152871 59 90
17 2 362301 76 110
18 -4 183167 91 138
19 3 277965 115 133
20 2 218946 76 96
21 2 244052 101 164
22 5 233328 92 102
23 -2 206161 75 99
24 -2 207176 56 114
25 -3 196553 41 99
26 2 143246 67 104
27 2 182192 77 138
28 2 194979 66 151
29 4 143756 105 120
30 4 275541 116 115
31 2 152299 62 98
32 2 193339 100 71
33 -4 130585 67 107
34 3 112611 46 73
35 3 148446 135 129
36 2 182079 124 118
37 -1 243060 58 104
38 -3 162765 68 107
39 1 225060 93 139
40 -3 133328 56 56
41 3 100750 83 93
42 3 132487 71 98
43 -3 317394 116 82
44 -4 184510 64 140
45 2 128423 32 120
46 -1 97839 25 66
47 3 172494 46 139
48 2 229242 63 119
49 5 351619 95 141
50 2 324598 113 133
51 -2 195838 111 98
52 3 199476 87 105
53 -2 92499 25 55
54 6 181633 47 73
55 -3 271856 109 86
56 3 95227 37 48
57 -2 118612 54 43
58 1 65475 16 46
59 2 121848 37 52
60 2 76302 29 68
61 -3 98104 55 47
62 -2 30989 5 41
63 1 31774 0 47
64 -4 150580 27 71
65 1 59382 29 24
> k <- length(x[1,])
> df <- as.data.frame(x)
> (mylm <- lm(df))
Call:
lm(formula = df)
Coefficients:
(Intercept) Time Comp LFB
-7.250e-01 4.431e-07 8.040e-03 7.678e-03
> (mysum <- summary(mylm))
Call:
lm(formula = df)
Residuals:
Min 1Q Median 3Q Max
-5.180 -2.692 1.128 1.974 5.706
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) -7.250e-01 1.159e+00 -0.625 0.534
Time 4.431e-07 6.653e-06 0.067 0.947
Comp 8.040e-03 1.667e-02 0.482 0.631
LFB 7.678e-03 1.493e-02 0.514 0.609
Residual standard error: 2.826 on 61 degrees of freedom
Multiple R-squared: 0.02927, Adjusted R-squared: -0.01847
F-statistic: 0.6131 on 3 and 61 DF, p-value: 0.6091
> 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.6109066 0.7781867 0.3890934
[2,] 0.5223628 0.9552743 0.4776372
[3,] 0.4450432 0.8900864 0.5549568
[4,] 0.5052377 0.9895247 0.4947623
[5,] 0.4471999 0.8943997 0.5528001
[6,] 0.5171830 0.9656340 0.4828170
[7,] 0.5645354 0.8709292 0.4354646
[8,] 0.7001891 0.5996218 0.2998109
[9,] 0.8132934 0.3734133 0.1867066
[10,] 0.7592205 0.4815590 0.2407795
[11,] 0.6893798 0.6212403 0.3106202
[12,] 0.7604323 0.4791353 0.2395677
[13,] 0.7342262 0.5315476 0.2657738
[14,] 0.6716859 0.6566282 0.3283141
[15,] 0.6341062 0.7317876 0.3658938
[16,] 0.6975036 0.6049927 0.3024964
[17,] 0.6934227 0.6131546 0.3065773
[18,] 0.6695132 0.6609736 0.3304868
[19,] 0.6763415 0.6473170 0.3236585
[20,] 0.6327521 0.7344957 0.3672479
[21,] 0.5925578 0.8148845 0.4074422
[22,] 0.5512907 0.8974187 0.4487093
[23,] 0.5423921 0.9152159 0.4576079
[24,] 0.5260521 0.9478958 0.4739479
[25,] 0.4684909 0.9369817 0.5315091
[26,] 0.4158022 0.8316044 0.5841978
[27,] 0.5443521 0.9112958 0.4556479
[28,] 0.5372428 0.9255143 0.4627572
[29,] 0.4793729 0.9587458 0.5206271
[30,] 0.4203438 0.8406876 0.5796562
[31,] 0.3733140 0.7466279 0.6266860
[32,] 0.4193617 0.8387234 0.5806383
[33,] 0.3473507 0.6947013 0.6526493
[34,] 0.3611565 0.7223129 0.6388435
[35,] 0.3650707 0.7301414 0.6349293
[36,] 0.3857254 0.7714507 0.6142746
[37,] 0.4723894 0.9447788 0.5276106
[38,] 0.6525502 0.6948996 0.3474498
[39,] 0.5951378 0.8097243 0.4048622
[40,] 0.5398943 0.9202114 0.4601057
[41,] 0.4829790 0.9659580 0.5170210
[42,] 0.4101279 0.8202559 0.5898721
[43,] 0.3784160 0.7568320 0.6215840
[44,] 0.2962588 0.5925176 0.7037412
[45,] 0.2611591 0.5223182 0.7388409
[46,] 0.2360056 0.4720112 0.7639944
[47,] 0.2012081 0.4024161 0.7987919
[48,] 0.5765644 0.8468711 0.4234356
[49,] 0.4799714 0.9599428 0.5200286
[50,] 0.5333000 0.9333999 0.4667000
[51,] 0.4026523 0.8053046 0.5973477
[52,] 0.2719352 0.5438704 0.7280648
> postscript(file="/var/wessaorg/rcomp/tmp/1j1wx1324590308.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/2bujn1324590308.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/335en1324590308.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/4bcjg1324590308.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/5g2gt1324590308.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 = 65
Frequency = 1
1 2 3 4 5 6 7
1.2746541 2.9863898 -5.1795859 3.6212269 2.8220613 -2.0803901 0.1000161
8 9 10 11 12 13 14
2.1855032 -2.1396511 2.7821886 1.9381553 1.1284421 -3.0859628 -3.9607759
15 16 17 18 19 20 21
-4.5207237 1.4918810 1.1088421 -5.1473660 1.6560601 1.2798558 0.5456277
22 23 24 25 26 27 28
4.0987760 -2.7294732 -2.6923345 -3.4518581 1.3243342 0.9656251 0.9485842
29 30 31 32 33 34 35
2.8957424 2.7872985 1.4065905 1.2901943 -4.6930897 2.7447658 1.5833640
36 37 38 39 40 41 42
0.7413585 -1.6475346 -3.7153888 -0.1896869 -3.2142871 2.2989836 2.3430099
43 44 45 46 47 48 49
-3.9778723 -4.9462389 1.4894518 -1.0261039 2.2114823 1.2032184 3.7227982
50 51 52 53 54 55 56
0.6514766 -3.0066587 2.1309416 -1.9392795 5.7061418 -3.9321267 3.0167786
57 58 59 60 61 62 63
-2.0918723 1.2141565 1.9742706 1.9359235 -3.1215371 -1.6437332 1.3500506
64 65
-4.1039437 1.2812535
> postscript(file="/var/wessaorg/rcomp/tmp/60i0a1324590308.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 = 65
Frequency = 1
lag(myerror, k = 1) myerror
0 1.2746541 NA
1 2.9863898 1.2746541
2 -5.1795859 2.9863898
3 3.6212269 -5.1795859
4 2.8220613 3.6212269
5 -2.0803901 2.8220613
6 0.1000161 -2.0803901
7 2.1855032 0.1000161
8 -2.1396511 2.1855032
9 2.7821886 -2.1396511
10 1.9381553 2.7821886
11 1.1284421 1.9381553
12 -3.0859628 1.1284421
13 -3.9607759 -3.0859628
14 -4.5207237 -3.9607759
15 1.4918810 -4.5207237
16 1.1088421 1.4918810
17 -5.1473660 1.1088421
18 1.6560601 -5.1473660
19 1.2798558 1.6560601
20 0.5456277 1.2798558
21 4.0987760 0.5456277
22 -2.7294732 4.0987760
23 -2.6923345 -2.7294732
24 -3.4518581 -2.6923345
25 1.3243342 -3.4518581
26 0.9656251 1.3243342
27 0.9485842 0.9656251
28 2.8957424 0.9485842
29 2.7872985 2.8957424
30 1.4065905 2.7872985
31 1.2901943 1.4065905
32 -4.6930897 1.2901943
33 2.7447658 -4.6930897
34 1.5833640 2.7447658
35 0.7413585 1.5833640
36 -1.6475346 0.7413585
37 -3.7153888 -1.6475346
38 -0.1896869 -3.7153888
39 -3.2142871 -0.1896869
40 2.2989836 -3.2142871
41 2.3430099 2.2989836
42 -3.9778723 2.3430099
43 -4.9462389 -3.9778723
44 1.4894518 -4.9462389
45 -1.0261039 1.4894518
46 2.2114823 -1.0261039
47 1.2032184 2.2114823
48 3.7227982 1.2032184
49 0.6514766 3.7227982
50 -3.0066587 0.6514766
51 2.1309416 -3.0066587
52 -1.9392795 2.1309416
53 5.7061418 -1.9392795
54 -3.9321267 5.7061418
55 3.0167786 -3.9321267
56 -2.0918723 3.0167786
57 1.2141565 -2.0918723
58 1.9742706 1.2141565
59 1.9359235 1.9742706
60 -3.1215371 1.9359235
61 -1.6437332 -3.1215371
62 1.3500506 -1.6437332
63 -4.1039437 1.3500506
64 1.2812535 -4.1039437
65 NA 1.2812535
> dum1 <- dum[2:length(myerror),]
> dum1
lag(myerror, k = 1) myerror
[1,] 2.9863898 1.2746541
[2,] -5.1795859 2.9863898
[3,] 3.6212269 -5.1795859
[4,] 2.8220613 3.6212269
[5,] -2.0803901 2.8220613
[6,] 0.1000161 -2.0803901
[7,] 2.1855032 0.1000161
[8,] -2.1396511 2.1855032
[9,] 2.7821886 -2.1396511
[10,] 1.9381553 2.7821886
[11,] 1.1284421 1.9381553
[12,] -3.0859628 1.1284421
[13,] -3.9607759 -3.0859628
[14,] -4.5207237 -3.9607759
[15,] 1.4918810 -4.5207237
[16,] 1.1088421 1.4918810
[17,] -5.1473660 1.1088421
[18,] 1.6560601 -5.1473660
[19,] 1.2798558 1.6560601
[20,] 0.5456277 1.2798558
[21,] 4.0987760 0.5456277
[22,] -2.7294732 4.0987760
[23,] -2.6923345 -2.7294732
[24,] -3.4518581 -2.6923345
[25,] 1.3243342 -3.4518581
[26,] 0.9656251 1.3243342
[27,] 0.9485842 0.9656251
[28,] 2.8957424 0.9485842
[29,] 2.7872985 2.8957424
[30,] 1.4065905 2.7872985
[31,] 1.2901943 1.4065905
[32,] -4.6930897 1.2901943
[33,] 2.7447658 -4.6930897
[34,] 1.5833640 2.7447658
[35,] 0.7413585 1.5833640
[36,] -1.6475346 0.7413585
[37,] -3.7153888 -1.6475346
[38,] -0.1896869 -3.7153888
[39,] -3.2142871 -0.1896869
[40,] 2.2989836 -3.2142871
[41,] 2.3430099 2.2989836
[42,] -3.9778723 2.3430099
[43,] -4.9462389 -3.9778723
[44,] 1.4894518 -4.9462389
[45,] -1.0261039 1.4894518
[46,] 2.2114823 -1.0261039
[47,] 1.2032184 2.2114823
[48,] 3.7227982 1.2032184
[49,] 0.6514766 3.7227982
[50,] -3.0066587 0.6514766
[51,] 2.1309416 -3.0066587
[52,] -1.9392795 2.1309416
[53,] 5.7061418 -1.9392795
[54,] -3.9321267 5.7061418
[55,] 3.0167786 -3.9321267
[56,] -2.0918723 3.0167786
[57,] 1.2141565 -2.0918723
[58,] 1.9742706 1.2141565
[59,] 1.9359235 1.9742706
[60,] -3.1215371 1.9359235
[61,] -1.6437332 -3.1215371
[62,] 1.3500506 -1.6437332
[63,] -4.1039437 1.3500506
[64,] 1.2812535 -4.1039437
> z <- as.data.frame(dum1)
> z
lag(myerror, k = 1) myerror
1 2.9863898 1.2746541
2 -5.1795859 2.9863898
3 3.6212269 -5.1795859
4 2.8220613 3.6212269
5 -2.0803901 2.8220613
6 0.1000161 -2.0803901
7 2.1855032 0.1000161
8 -2.1396511 2.1855032
9 2.7821886 -2.1396511
10 1.9381553 2.7821886
11 1.1284421 1.9381553
12 -3.0859628 1.1284421
13 -3.9607759 -3.0859628
14 -4.5207237 -3.9607759
15 1.4918810 -4.5207237
16 1.1088421 1.4918810
17 -5.1473660 1.1088421
18 1.6560601 -5.1473660
19 1.2798558 1.6560601
20 0.5456277 1.2798558
21 4.0987760 0.5456277
22 -2.7294732 4.0987760
23 -2.6923345 -2.7294732
24 -3.4518581 -2.6923345
25 1.3243342 -3.4518581
26 0.9656251 1.3243342
27 0.9485842 0.9656251
28 2.8957424 0.9485842
29 2.7872985 2.8957424
30 1.4065905 2.7872985
31 1.2901943 1.4065905
32 -4.6930897 1.2901943
33 2.7447658 -4.6930897
34 1.5833640 2.7447658
35 0.7413585 1.5833640
36 -1.6475346 0.7413585
37 -3.7153888 -1.6475346
38 -0.1896869 -3.7153888
39 -3.2142871 -0.1896869
40 2.2989836 -3.2142871
41 2.3430099 2.2989836
42 -3.9778723 2.3430099
43 -4.9462389 -3.9778723
44 1.4894518 -4.9462389
45 -1.0261039 1.4894518
46 2.2114823 -1.0261039
47 1.2032184 2.2114823
48 3.7227982 1.2032184
49 0.6514766 3.7227982
50 -3.0066587 0.6514766
51 2.1309416 -3.0066587
52 -1.9392795 2.1309416
53 5.7061418 -1.9392795
54 -3.9321267 5.7061418
55 3.0167786 -3.9321267
56 -2.0918723 3.0167786
57 1.2141565 -2.0918723
58 1.9742706 1.2141565
59 1.9359235 1.9742706
60 -3.1215371 1.9359235
61 -1.6437332 -3.1215371
62 1.3500506 -1.6437332
63 -4.1039437 1.3500506
64 1.2812535 -4.1039437
> 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/7b68u1324590308.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/8n6661324590308.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/9733j1324590308.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/10xn611324590308.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/11slo41324590308.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/12x8le1324590308.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/13b17b1324590308.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/14j7lv1324590308.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/15duq01324590308.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/160x781324590308.tab")
+ }
>
> try(system("convert tmp/1j1wx1324590308.ps tmp/1j1wx1324590308.png",intern=TRUE))
character(0)
> try(system("convert tmp/2bujn1324590308.ps tmp/2bujn1324590308.png",intern=TRUE))
character(0)
> try(system("convert tmp/335en1324590308.ps tmp/335en1324590308.png",intern=TRUE))
character(0)
> try(system("convert tmp/4bcjg1324590308.ps tmp/4bcjg1324590308.png",intern=TRUE))
character(0)
> try(system("convert tmp/5g2gt1324590308.ps tmp/5g2gt1324590308.png",intern=TRUE))
character(0)
> try(system("convert tmp/60i0a1324590308.ps tmp/60i0a1324590308.png",intern=TRUE))
character(0)
> try(system("convert tmp/7b68u1324590308.ps tmp/7b68u1324590308.png",intern=TRUE))
character(0)
> try(system("convert tmp/8n6661324590308.ps tmp/8n6661324590308.png",intern=TRUE))
character(0)
> try(system("convert tmp/9733j1324590308.ps tmp/9733j1324590308.png",intern=TRUE))
character(0)
> try(system("convert tmp/10xn611324590308.ps tmp/10xn611324590308.png",intern=TRUE))
character(0)
>
>
> proc.time()
user system elapsed
3.450 0.676 4.139