R version 2.9.0 (2009-04-17)
Copyright (C) 2009 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.
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(2756.76
+ ,10001.60
+ ,2849.27
+ ,10411.75
+ ,2921.44
+ ,10673.38
+ ,2981.85
+ ,10539.51
+ ,3080.58
+ ,10723.78
+ ,3106.22
+ ,10682.06
+ ,3119.31
+ ,10283.19
+ ,3061.26
+ ,10377.18
+ ,3097.31
+ ,10486.64
+ ,3161.69
+ ,10545.38
+ ,3257.16
+ ,10554.27
+ ,3277.01
+ ,10532.54
+ ,3295.32
+ ,10324.31
+ ,3363.99
+ ,10695.25
+ ,3494.17
+ ,10827.81
+ ,3667.03
+ ,10872.48
+ ,3813.06
+ ,10971.19
+ ,3917.96
+ ,11145.65
+ ,3895.51
+ ,11234.68
+ ,3801.06
+ ,11333.88
+ ,3570.12
+ ,10997.97
+ ,3701.61
+ ,11036.89
+ ,3862.27
+ ,11257.35
+ ,3970.10
+ ,11533.59
+ ,4138.52
+ ,11963.12
+ ,4199.75
+ ,12185.15
+ ,4290.89
+ ,12377.62
+ ,4443.91
+ ,12512.89
+ ,4502.64
+ ,12631.48
+ ,4356.98
+ ,12268.53
+ ,4591.27
+ ,12754.80
+ ,4696.96
+ ,13407.75
+ ,4621.40
+ ,13480.21
+ ,4562.84
+ ,13673.28
+ ,4202.52
+ ,13239.71
+ ,4296.49
+ ,13557.69
+ ,4435.23
+ ,13901.28
+ ,4105.18
+ ,13200.58
+ ,4116.68
+ ,13406.97
+ ,3844.49
+ ,12538.12
+ ,3720.98
+ ,12419.57
+ ,3674.40
+ ,12193.88
+ ,3857.62
+ ,12656.63
+ ,3801.06
+ ,12812.48
+ ,3504.37
+ ,12056.67
+ ,3032.60
+ ,11322.38
+ ,3047.03
+ ,11530.75
+ ,2962.34
+ ,11114.08
+ ,2197.82
+ ,9181.73
+ ,2014.45
+ ,8614.55
+ ,1862.83
+ ,8595.56
+ ,1905.41
+ ,8396.20
+ ,1810.99
+ ,7690.50
+ ,1670.07
+ ,7235.47
+ ,1864.44
+ ,7992.12
+ ,2052.02
+ ,8398.37
+ ,2029.60
+ ,8593.01
+ ,2070.83
+ ,8679.75
+ ,2293.41
+ ,9374.63
+ ,2443.27
+ ,9634.97)
+ ,dim=c(2
+ ,60)
+ ,dimnames=list(c('Bel20'
+ ,'Dow
')
+ ,1:60))
> y <- array(NA,dim=c(2,60),dimnames=list(c('Bel20','Dow
'),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 = 'Include Monthly Dummies'
> par1 = '2'
> #'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
Dow\r Bel20 M1 M2 M3 M4 M5 M6 M7 M8 M9 M10 M11 t
1 10001.60 2756.76 1 0 0 0 0 0 0 0 0 0 0 1
2 10411.75 2849.27 0 1 0 0 0 0 0 0 0 0 0 2
3 10673.38 2921.44 0 0 1 0 0 0 0 0 0 0 0 3
4 10539.51 2981.85 0 0 0 1 0 0 0 0 0 0 0 4
5 10723.78 3080.58 0 0 0 0 1 0 0 0 0 0 0 5
6 10682.06 3106.22 0 0 0 0 0 1 0 0 0 0 0 6
7 10283.19 3119.31 0 0 0 0 0 0 1 0 0 0 0 7
8 10377.18 3061.26 0 0 0 0 0 0 0 1 0 0 0 8
9 10486.64 3097.31 0 0 0 0 0 0 0 0 1 0 0 9
10 10545.38 3161.69 0 0 0 0 0 0 0 0 0 1 0 10
11 10554.27 3257.16 0 0 0 0 0 0 0 0 0 0 1 11
12 10532.54 3277.01 0 0 0 0 0 0 0 0 0 0 0 12
13 10324.31 3295.32 1 0 0 0 0 0 0 0 0 0 0 13
14 10695.25 3363.99 0 1 0 0 0 0 0 0 0 0 0 14
15 10827.81 3494.17 0 0 1 0 0 0 0 0 0 0 0 15
16 10872.48 3667.03 0 0 0 1 0 0 0 0 0 0 0 16
17 10971.19 3813.06 0 0 0 0 1 0 0 0 0 0 0 17
18 11145.65 3917.96 0 0 0 0 0 1 0 0 0 0 0 18
19 11234.68 3895.51 0 0 0 0 0 0 1 0 0 0 0 19
20 11333.88 3801.06 0 0 0 0 0 0 0 1 0 0 0 20
21 10997.97 3570.12 0 0 0 0 0 0 0 0 1 0 0 21
22 11036.89 3701.61 0 0 0 0 0 0 0 0 0 1 0 22
23 11257.35 3862.27 0 0 0 0 0 0 0 0 0 0 1 23
24 11533.59 3970.10 0 0 0 0 0 0 0 0 0 0 0 24
25 11963.12 4138.52 1 0 0 0 0 0 0 0 0 0 0 25
26 12185.15 4199.75 0 1 0 0 0 0 0 0 0 0 0 26
27 12377.62 4290.89 0 0 1 0 0 0 0 0 0 0 0 27
28 12512.89 4443.91 0 0 0 1 0 0 0 0 0 0 0 28
29 12631.48 4502.64 0 0 0 0 1 0 0 0 0 0 0 29
30 12268.53 4356.98 0 0 0 0 0 1 0 0 0 0 0 30
31 12754.80 4591.27 0 0 0 0 0 0 1 0 0 0 0 31
32 13407.75 4696.96 0 0 0 0 0 0 0 1 0 0 0 32
33 13480.21 4621.40 0 0 0 0 0 0 0 0 1 0 0 33
34 13673.28 4562.84 0 0 0 0 0 0 0 0 0 1 0 34
35 13239.71 4202.52 0 0 0 0 0 0 0 0 0 0 1 35
36 13557.69 4296.49 0 0 0 0 0 0 0 0 0 0 0 36
37 13901.28 4435.23 1 0 0 0 0 0 0 0 0 0 0 37
38 13200.58 4105.18 0 1 0 0 0 0 0 0 0 0 0 38
39 13406.97 4116.68 0 0 1 0 0 0 0 0 0 0 0 39
40 12538.12 3844.49 0 0 0 1 0 0 0 0 0 0 0 40
41 12419.57 3720.98 0 0 0 0 1 0 0 0 0 0 0 41
42 12193.88 3674.40 0 0 0 0 0 1 0 0 0 0 0 42
43 12656.63 3857.62 0 0 0 0 0 0 1 0 0 0 0 43
44 12812.48 3801.06 0 0 0 0 0 0 0 1 0 0 0 44
45 12056.67 3504.37 0 0 0 0 0 0 0 0 1 0 0 45
46 11322.38 3032.60 0 0 0 0 0 0 0 0 0 1 0 46
47 11530.75 3047.03 0 0 0 0 0 0 0 0 0 0 1 47
48 11114.08 2962.34 0 0 0 0 0 0 0 0 0 0 0 48
49 9181.73 2197.82 1 0 0 0 0 0 0 0 0 0 0 49
50 8614.55 2014.45 0 1 0 0 0 0 0 0 0 0 0 50
51 8595.56 1862.83 0 0 1 0 0 0 0 0 0 0 0 51
52 8396.20 1905.41 0 0 0 1 0 0 0 0 0 0 0 52
53 7690.50 1810.99 0 0 0 0 1 0 0 0 0 0 0 53
54 7235.47 1670.07 0 0 0 0 0 1 0 0 0 0 0 54
55 7992.12 1864.44 0 0 0 0 0 0 1 0 0 0 0 55
56 8398.37 2052.02 0 0 0 0 0 0 0 1 0 0 0 56
57 8593.01 2029.60 0 0 0 0 0 0 0 0 1 0 0 57
58 8679.75 2070.83 0 0 0 0 0 0 0 0 0 1 0 58
59 9374.63 2293.41 0 0 0 0 0 0 0 0 0 0 1 59
60 9634.97 2443.27 0 0 0 0 0 0 0 0 0 0 0 60
> k <- length(x[1,])
> df <- as.data.frame(x)
> (mylm <- lm(df))
Call:
lm(formula = df)
Coefficients:
(Intercept) Bel20 M1 M2 M3 M4
4315.150 1.918 -12.011 33.934 117.190 -160.063
M5 M6 M7 M8 M9 M10
-290.144 -407.337 -372.008 -135.387 -65.007 -36.622
M11 t
39.511 12.727
> (mysum <- summary(mylm))
Call:
lm(formula = df)
Residuals:
Min 1Q Median 3Q Max
-797.3 -423.4 -129.5 508.2 782.8
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 4315.15037 468.49562 9.211 5.2e-12 ***
Bel20 1.91787 0.09126 21.016 < 2e-16 ***
M1 -12.01144 355.21187 -0.034 0.97317
M2 33.93358 354.91003 0.096 0.92424
M3 117.19030 354.16622 0.331 0.74223
M4 -160.06261 353.52353 -0.453 0.65285
M5 -290.14379 353.02662 -0.822 0.41539
M6 -407.33681 352.77229 -1.155 0.25419
M7 -372.00804 352.23065 -1.056 0.29642
M8 -135.38739 352.00662 -0.385 0.70229
M9 -65.00656 351.89407 -0.185 0.85425
M10 -36.62201 351.86278 -0.104 0.91756
M11 39.51114 351.68555 0.112 0.91104
t 12.72665 4.51960 2.816 0.00714 **
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 555.9 on 46 degrees of freedom
Multiple R-squared: 0.911, Adjusted R-squared: 0.8858
F-statistic: 36.2 on 13 and 46 DF, p-value: < 2.2e-16
> 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.003299534 0.0065990679 0.9967004660
[2,] 0.001055753 0.0021115051 0.9989442475
[3,] 0.017977261 0.0359545225 0.9820227387
[4,] 0.023534049 0.0470680973 0.9764659513
[5,] 0.020726487 0.0414529742 0.9792735129
[6,] 0.010480340 0.0209606795 0.9895196603
[7,] 0.006008196 0.0120163910 0.9939918045
[8,] 0.006800462 0.0136009248 0.9931995376
[9,] 0.027503867 0.0550077332 0.9724961334
[10,] 0.023075273 0.0461505456 0.9769247272
[11,] 0.020331270 0.0406625402 0.9796687299
[12,] 0.021835730 0.0436714606 0.9781642697
[13,] 0.028273388 0.0565467760 0.9717266120
[14,] 0.037742790 0.0754855799 0.9622572100
[15,] 0.089875229 0.1797504570 0.9101247715
[16,] 0.170643425 0.3412868491 0.8293565755
[17,] 0.296186993 0.5923739851 0.7038130075
[18,] 0.672067978 0.6558640445 0.3279320222
[19,] 0.969116333 0.0617673346 0.0308836673
[20,] 0.998180851 0.0036382985 0.0018191493
[21,] 0.998284821 0.0034303573 0.0017151787
[22,] 0.996271272 0.0074574562 0.0037287281
[23,] 0.996280071 0.0074398578 0.0037199289
[24,] 0.999527098 0.0009458037 0.0004729019
[25,] 0.998059920 0.0038801594 0.0019400797
[26,] 0.993999399 0.0120012013 0.0060006006
[27,] 0.978108554 0.0437828926 0.0218914463
> postscript(file="/var/www/html/rcomp/tmp/146o01259619073.ps",horizontal=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/rcomp/tmp/2rwzi1259619073.ps",horizontal=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/rcomp/tmp/3z6ek1259619073.ps",horizontal=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/rcomp/tmp/4yctd1259619073.ps",horizontal=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/rcomp/tmp/5gsw21259619073.ps",horizontal=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 7
398.63678 572.69328 599.92749 614.72544 726.99902 740.57130 268.54101
8 9 10 11 12 13 14
224.51587 181.72930 75.88586 -187.18265 -220.19780 -464.25914 -283.69070
15 16 17 18 19 20 21
-496.78193 -519.10806 -583.10958 -505.36740 -421.33671 -390.34151 -366.44690
22 23 24 25 26 27 28
-620.81836 -797.34259 -701.12164 -595.31392 -549.38655 -627.69426 -521.36993
29 30 31 32 33 34 35
-398.06170 -377.18889 -388.31124 -187.40785 -53.14133 211.12773 379.74360
36 37 38 39 40 41 42
544.28618 621.07617 494.69634 583.04752 500.74786 736.42809 704.53869
43 44 45 46 47 48 49
767.84177 782.81901 512.91334 642.30402 734.13941 506.67802 39.86011
50 51 52 53 54 55 56
-234.31237 -58.49882 -74.99531 -482.25582 -562.55370 -226.73483 -429.58552
57 58 59 60
-275.05442 -308.49926 -129.35778 -129.64476
> postscript(file="/var/www/html/rcomp/tmp/6g6yq1259619073.ps",horizontal=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 398.63678 NA
1 572.69328 398.63678
2 599.92749 572.69328
3 614.72544 599.92749
4 726.99902 614.72544
5 740.57130 726.99902
6 268.54101 740.57130
7 224.51587 268.54101
8 181.72930 224.51587
9 75.88586 181.72930
10 -187.18265 75.88586
11 -220.19780 -187.18265
12 -464.25914 -220.19780
13 -283.69070 -464.25914
14 -496.78193 -283.69070
15 -519.10806 -496.78193
16 -583.10958 -519.10806
17 -505.36740 -583.10958
18 -421.33671 -505.36740
19 -390.34151 -421.33671
20 -366.44690 -390.34151
21 -620.81836 -366.44690
22 -797.34259 -620.81836
23 -701.12164 -797.34259
24 -595.31392 -701.12164
25 -549.38655 -595.31392
26 -627.69426 -549.38655
27 -521.36993 -627.69426
28 -398.06170 -521.36993
29 -377.18889 -398.06170
30 -388.31124 -377.18889
31 -187.40785 -388.31124
32 -53.14133 -187.40785
33 211.12773 -53.14133
34 379.74360 211.12773
35 544.28618 379.74360
36 621.07617 544.28618
37 494.69634 621.07617
38 583.04752 494.69634
39 500.74786 583.04752
40 736.42809 500.74786
41 704.53869 736.42809
42 767.84177 704.53869
43 782.81901 767.84177
44 512.91334 782.81901
45 642.30402 512.91334
46 734.13941 642.30402
47 506.67802 734.13941
48 39.86011 506.67802
49 -234.31237 39.86011
50 -58.49882 -234.31237
51 -74.99531 -58.49882
52 -482.25582 -74.99531
53 -562.55370 -482.25582
54 -226.73483 -562.55370
55 -429.58552 -226.73483
56 -275.05442 -429.58552
57 -308.49926 -275.05442
58 -129.35778 -308.49926
59 -129.64476 -129.35778
60 NA -129.64476
> dum1 <- dum[2:length(myerror),]
> dum1
lag(myerror, k = 1) myerror
[1,] 572.69328 398.63678
[2,] 599.92749 572.69328
[3,] 614.72544 599.92749
[4,] 726.99902 614.72544
[5,] 740.57130 726.99902
[6,] 268.54101 740.57130
[7,] 224.51587 268.54101
[8,] 181.72930 224.51587
[9,] 75.88586 181.72930
[10,] -187.18265 75.88586
[11,] -220.19780 -187.18265
[12,] -464.25914 -220.19780
[13,] -283.69070 -464.25914
[14,] -496.78193 -283.69070
[15,] -519.10806 -496.78193
[16,] -583.10958 -519.10806
[17,] -505.36740 -583.10958
[18,] -421.33671 -505.36740
[19,] -390.34151 -421.33671
[20,] -366.44690 -390.34151
[21,] -620.81836 -366.44690
[22,] -797.34259 -620.81836
[23,] -701.12164 -797.34259
[24,] -595.31392 -701.12164
[25,] -549.38655 -595.31392
[26,] -627.69426 -549.38655
[27,] -521.36993 -627.69426
[28,] -398.06170 -521.36993
[29,] -377.18889 -398.06170
[30,] -388.31124 -377.18889
[31,] -187.40785 -388.31124
[32,] -53.14133 -187.40785
[33,] 211.12773 -53.14133
[34,] 379.74360 211.12773
[35,] 544.28618 379.74360
[36,] 621.07617 544.28618
[37,] 494.69634 621.07617
[38,] 583.04752 494.69634
[39,] 500.74786 583.04752
[40,] 736.42809 500.74786
[41,] 704.53869 736.42809
[42,] 767.84177 704.53869
[43,] 782.81901 767.84177
[44,] 512.91334 782.81901
[45,] 642.30402 512.91334
[46,] 734.13941 642.30402
[47,] 506.67802 734.13941
[48,] 39.86011 506.67802
[49,] -234.31237 39.86011
[50,] -58.49882 -234.31237
[51,] -74.99531 -58.49882
[52,] -482.25582 -74.99531
[53,] -562.55370 -482.25582
[54,] -226.73483 -562.55370
[55,] -429.58552 -226.73483
[56,] -275.05442 -429.58552
[57,] -308.49926 -275.05442
[58,] -129.35778 -308.49926
[59,] -129.64476 -129.35778
> z <- as.data.frame(dum1)
> z
lag(myerror, k = 1) myerror
1 572.69328 398.63678
2 599.92749 572.69328
3 614.72544 599.92749
4 726.99902 614.72544
5 740.57130 726.99902
6 268.54101 740.57130
7 224.51587 268.54101
8 181.72930 224.51587
9 75.88586 181.72930
10 -187.18265 75.88586
11 -220.19780 -187.18265
12 -464.25914 -220.19780
13 -283.69070 -464.25914
14 -496.78193 -283.69070
15 -519.10806 -496.78193
16 -583.10958 -519.10806
17 -505.36740 -583.10958
18 -421.33671 -505.36740
19 -390.34151 -421.33671
20 -366.44690 -390.34151
21 -620.81836 -366.44690
22 -797.34259 -620.81836
23 -701.12164 -797.34259
24 -595.31392 -701.12164
25 -549.38655 -595.31392
26 -627.69426 -549.38655
27 -521.36993 -627.69426
28 -398.06170 -521.36993
29 -377.18889 -398.06170
30 -388.31124 -377.18889
31 -187.40785 -388.31124
32 -53.14133 -187.40785
33 211.12773 -53.14133
34 379.74360 211.12773
35 544.28618 379.74360
36 621.07617 544.28618
37 494.69634 621.07617
38 583.04752 494.69634
39 500.74786 583.04752
40 736.42809 500.74786
41 704.53869 736.42809
42 767.84177 704.53869
43 782.81901 767.84177
44 512.91334 782.81901
45 642.30402 512.91334
46 734.13941 642.30402
47 506.67802 734.13941
48 39.86011 506.67802
49 -234.31237 39.86011
50 -58.49882 -234.31237
51 -74.99531 -58.49882
52 -482.25582 -74.99531
53 -562.55370 -482.25582
54 -226.73483 -562.55370
55 -429.58552 -226.73483
56 -275.05442 -429.58552
57 -308.49926 -275.05442
58 -129.35778 -308.49926
59 -129.64476 -129.35778
> 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/rcomp/tmp/7bkbr1259619073.ps",horizontal=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/rcomp/tmp/8eua71259619073.ps",horizontal=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/rcomp/tmp/9ksuo1259619073.ps",horizontal=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/rcomp/tmp/10zeeu1259619073.ps",horizontal=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/rcomp/createtable file can be downloaded at http://www.wessa.net/cretab
> load(file="/var/www/html/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/rcomp/tmp/11rli51259619073.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/rcomp/tmp/12rbtm1259619073.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/rcomp/tmp/13h2pf1259619073.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/rcomp/tmp/14onrh1259619073.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/rcomp/tmp/156fnc1259619073.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/rcomp/tmp/16ljty1259619073.tab")
+ }
>
> system("convert tmp/146o01259619073.ps tmp/146o01259619073.png")
> system("convert tmp/2rwzi1259619073.ps tmp/2rwzi1259619073.png")
> system("convert tmp/3z6ek1259619073.ps tmp/3z6ek1259619073.png")
> system("convert tmp/4yctd1259619073.ps tmp/4yctd1259619073.png")
> system("convert tmp/5gsw21259619073.ps tmp/5gsw21259619073.png")
> system("convert tmp/6g6yq1259619073.ps tmp/6g6yq1259619073.png")
> system("convert tmp/7bkbr1259619073.ps tmp/7bkbr1259619073.png")
> system("convert tmp/8eua71259619073.ps tmp/8eua71259619073.png")
> system("convert tmp/9ksuo1259619073.ps tmp/9ksuo1259619073.png")
> system("convert tmp/10zeeu1259619073.ps tmp/10zeeu1259619073.png")
>
>
> proc.time()
user system elapsed
2.471 1.568 5.728