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(3016.70
+ ,2756.76
+ ,3052.40
+ ,2849.27
+ ,3099.60
+ ,2921.44
+ ,3103.30
+ ,2981.85
+ ,3119.80
+ ,3080.58
+ ,3093.70
+ ,3106.22
+ ,3164.90
+ ,3119.31
+ ,3311.50
+ ,3061.26
+ ,3410.60
+ ,3097.31
+ ,3392.60
+ ,3161.69
+ ,3338.20
+ ,3257.16
+ ,3285.10
+ ,3277.01
+ ,3294.80
+ ,3295.32
+ ,3611.20
+ ,3363.99
+ ,3611.30
+ ,3494.17
+ ,3521.00
+ ,3667.03
+ ,3519.30
+ ,3813.06
+ ,3438.30
+ ,3917.96
+ ,3534.90
+ ,3895.51
+ ,3705.80
+ ,3801.06
+ ,3807.60
+ ,3570.12
+ ,3663.00
+ ,3701.61
+ ,3604.50
+ ,3862.27
+ ,3563.80
+ ,3970.10
+ ,3511.40
+ ,4138.52
+ ,3546.50
+ ,4199.75
+ ,3525.40
+ ,4290.89
+ ,3529.90
+ ,4443.91
+ ,3591.60
+ ,4502.64
+ ,3668.30
+ ,4356.98
+ ,3728.80
+ ,4591.27
+ ,3853.60
+ ,4696.96
+ ,3897.70
+ ,4621.40
+ ,3640.70
+ ,4562.84
+ ,3495.50
+ ,4202.52
+ ,3495.10
+ ,4296.49
+ ,3268.00
+ ,4435.23
+ ,3479.10
+ ,4105.18
+ ,3417.80
+ ,4116.68
+ ,3521.30
+ ,3844.49
+ ,3487.10
+ ,3720.98
+ ,3529.90
+ ,3674.40
+ ,3544.30
+ ,3857.62
+ ,3710.80
+ ,3801.06
+ ,3641.90
+ ,3504.37
+ ,3447.10
+ ,3032.60
+ ,3386.80
+ ,3047.03
+ ,3438.50
+ ,2962.34
+ ,3364.30
+ ,2197.82
+ ,3462.70
+ ,2014.45
+ ,3291.90
+ ,1862.83
+ ,3550.00
+ ,1905.41
+ ,3611.00
+ ,1810.99
+ ,3708.60
+ ,1670.07
+ ,3771.10
+ ,1864.44
+ ,4042.70
+ ,2052.02
+ ,3988.40
+ ,2029.60
+ ,3851.20
+ ,2070.83
+ ,3876.70
+ ,2293.41)
+ ,dim=c(2
+ ,59)
+ ,dimnames=list(c('Zichtrekeningen'
+ ,'Bel20
')
+ ,1:59))
> y <- array(NA,dim=c(2,59),dimnames=list(c('Zichtrekeningen','Bel20
'),1:59))
> 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 = '1'
> #'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
Zichtrekeningen Bel20\r M1 M2 M3 M4 M5 M6 M7 M8 M9 M10 M11 t
1 3016.7 2756.76 1 0 0 0 0 0 0 0 0 0 0 1
2 3052.4 2849.27 0 1 0 0 0 0 0 0 0 0 0 2
3 3099.6 2921.44 0 0 1 0 0 0 0 0 0 0 0 3
4 3103.3 2981.85 0 0 0 1 0 0 0 0 0 0 0 4
5 3119.8 3080.58 0 0 0 0 1 0 0 0 0 0 0 5
6 3093.7 3106.22 0 0 0 0 0 1 0 0 0 0 0 6
7 3164.9 3119.31 0 0 0 0 0 0 1 0 0 0 0 7
8 3311.5 3061.26 0 0 0 0 0 0 0 1 0 0 0 8
9 3410.6 3097.31 0 0 0 0 0 0 0 0 1 0 0 9
10 3392.6 3161.69 0 0 0 0 0 0 0 0 0 1 0 10
11 3338.2 3257.16 0 0 0 0 0 0 0 0 0 0 1 11
12 3285.1 3277.01 0 0 0 0 0 0 0 0 0 0 0 12
13 3294.8 3295.32 1 0 0 0 0 0 0 0 0 0 0 13
14 3611.2 3363.99 0 1 0 0 0 0 0 0 0 0 0 14
15 3611.3 3494.17 0 0 1 0 0 0 0 0 0 0 0 15
16 3521.0 3667.03 0 0 0 1 0 0 0 0 0 0 0 16
17 3519.3 3813.06 0 0 0 0 1 0 0 0 0 0 0 17
18 3438.3 3917.96 0 0 0 0 0 1 0 0 0 0 0 18
19 3534.9 3895.51 0 0 0 0 0 0 1 0 0 0 0 19
20 3705.8 3801.06 0 0 0 0 0 0 0 1 0 0 0 20
21 3807.6 3570.12 0 0 0 0 0 0 0 0 1 0 0 21
22 3663.0 3701.61 0 0 0 0 0 0 0 0 0 1 0 22
23 3604.5 3862.27 0 0 0 0 0 0 0 0 0 0 1 23
24 3563.8 3970.10 0 0 0 0 0 0 0 0 0 0 0 24
25 3511.4 4138.52 1 0 0 0 0 0 0 0 0 0 0 25
26 3546.5 4199.75 0 1 0 0 0 0 0 0 0 0 0 26
27 3525.4 4290.89 0 0 1 0 0 0 0 0 0 0 0 27
28 3529.9 4443.91 0 0 0 1 0 0 0 0 0 0 0 28
29 3591.6 4502.64 0 0 0 0 1 0 0 0 0 0 0 29
30 3668.3 4356.98 0 0 0 0 0 1 0 0 0 0 0 30
31 3728.8 4591.27 0 0 0 0 0 0 1 0 0 0 0 31
32 3853.6 4696.96 0 0 0 0 0 0 0 1 0 0 0 32
33 3897.7 4621.40 0 0 0 0 0 0 0 0 1 0 0 33
34 3640.7 4562.84 0 0 0 0 0 0 0 0 0 1 0 34
35 3495.5 4202.52 0 0 0 0 0 0 0 0 0 0 1 35
36 3495.1 4296.49 0 0 0 0 0 0 0 0 0 0 0 36
37 3268.0 4435.23 1 0 0 0 0 0 0 0 0 0 0 37
38 3479.1 4105.18 0 1 0 0 0 0 0 0 0 0 0 38
39 3417.8 4116.68 0 0 1 0 0 0 0 0 0 0 0 39
40 3521.3 3844.49 0 0 0 1 0 0 0 0 0 0 0 40
41 3487.1 3720.98 0 0 0 0 1 0 0 0 0 0 0 41
42 3529.9 3674.40 0 0 0 0 0 1 0 0 0 0 0 42
43 3544.3 3857.62 0 0 0 0 0 0 1 0 0 0 0 43
44 3710.8 3801.06 0 0 0 0 0 0 0 1 0 0 0 44
45 3641.9 3504.37 0 0 0 0 0 0 0 0 1 0 0 45
46 3447.1 3032.60 0 0 0 0 0 0 0 0 0 1 0 46
47 3386.8 3047.03 0 0 0 0 0 0 0 0 0 0 1 47
48 3438.5 2962.34 0 0 0 0 0 0 0 0 0 0 0 48
49 3364.3 2197.82 1 0 0 0 0 0 0 0 0 0 0 49
50 3462.7 2014.45 0 1 0 0 0 0 0 0 0 0 0 50
51 3291.9 1862.83 0 0 1 0 0 0 0 0 0 0 0 51
52 3550.0 1905.41 0 0 0 1 0 0 0 0 0 0 0 52
53 3611.0 1810.99 0 0 0 0 1 0 0 0 0 0 0 53
54 3708.6 1670.07 0 0 0 0 0 1 0 0 0 0 0 54
55 3771.1 1864.44 0 0 0 0 0 0 1 0 0 0 0 55
56 4042.7 2052.02 0 0 0 0 0 0 0 1 0 0 0 56
57 3988.4 2029.60 0 0 0 0 0 0 0 0 1 0 0 57
58 3851.2 2070.83 0 0 0 0 0 0 0 0 0 1 0 58
59 3876.7 2293.41 0 0 0 0 0 0 0 0 0 0 1 59
> k <- length(x[1,])
> df <- as.data.frame(x)
> (mylm <- lm(df))
Call:
lm(formula = df)
Coefficients:
(Intercept) `Bel20\r` M1 M2 M3 M4
2943.37798 0.06832 -94.28827 40.54557 -11.21277 34.06367
M5 M6 M7 M8 M9 M10
45.07192 61.35797 105.68229 272.12898 296.06232 141.26650
M11 t
72.38897 8.48263
> (mysum <- summary(mylm))
Call:
lm(formula = df)
Residuals:
Min 1Q Median 3Q Max
-235.829 -113.117 7.117 118.769 313.167
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 2943.37798 133.93070 21.977 < 2e-16 ***
`Bel20\r` 0.06832 0.02523 2.708 0.00954 **
M1 -94.28827 103.14426 -0.914 0.36552
M2 40.54557 103.17658 0.393 0.69620
M3 -11.21277 103.02300 -0.109 0.91382
M4 34.06367 102.89552 0.331 0.74214
M5 45.07192 102.81785 0.438 0.66322
M6 61.35797 102.85101 0.597 0.55378
M7 105.68229 102.67721 1.029 0.30885
M8 272.12898 102.67144 2.650 0.01105 *
M9 296.06232 102.80740 2.880 0.00607 **
M10 141.26650 102.91634 1.373 0.17667
M11 72.38897 102.91603 0.703 0.48544
t 8.48263 1.25954 6.735 2.52e-08 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 153 on 45 degrees of freedom
Multiple R-squared: 0.6585, Adjusted R-squared: 0.5599
F-statistic: 6.675 on 13 and 45 DF, p-value: 7.135e-07
> 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.26656486 0.533129724 0.733435138
[2,] 0.14358406 0.287168119 0.856415941
[3,] 0.06981163 0.139623265 0.930188368
[4,] 0.03453023 0.069060452 0.965469774
[5,] 0.01866117 0.037322333 0.981338833
[6,] 0.02367147 0.047342943 0.976328529
[7,] 0.02066993 0.041339865 0.979330068
[8,] 0.01312167 0.026243333 0.986878334
[9,] 0.02092640 0.041852809 0.979073596
[10,] 0.06372208 0.127444150 0.936277925
[11,] 0.13851808 0.277036155 0.861481922
[12,] 0.11766501 0.235330014 0.882334993
[13,] 0.09711450 0.194229006 0.902885497
[14,] 0.08394432 0.167888645 0.916055677
[15,] 0.07547355 0.150947100 0.924526450
[16,] 0.06244536 0.124890711 0.937554645
[17,] 0.09439202 0.188784047 0.905607977
[18,] 0.20185506 0.403710119 0.798144941
[19,] 0.70764401 0.584711980 0.292355990
[20,] 0.96472057 0.070558865 0.035279433
[21,] 0.96831590 0.063368203 0.031684101
[22,] 0.95083300 0.098333991 0.049166996
[23,] 0.97415543 0.051689136 0.025844568
[24,] 0.99227258 0.015454839 0.007727419
[25,] 0.99598503 0.008029933 0.004014966
[26,] 0.99282471 0.014350583 0.007175292
> postscript(file="/var/www/html/rcomp/tmp/19hva1258742090.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/2ympw1258742090.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/3sfu31258742090.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/4687m1258742090.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/5emfa1258742090.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 = 59
Frequency = 1
1 2 3 4 5 6
-29.219301 -143.156225 -57.611309 -111.797698 -121.533992 -174.154439
7 8 9 10 11 12
-156.655715 -181.018954 -116.797925 7.116710 6.588924 16.039076
13 14 15 16 17 18
110.293749 278.685619 313.167184 157.298001 126.130083 13.194445
19 20 21 22 23 24
58.521328 60.945005 146.107288 138.836843 129.755155 145.594349
25 26 27 28 29 30
167.493228 55.093412 71.042262 11.328585 49.525166 111.408253
31 32 33 34 35 36
103.094182 45.743921 62.590354 -44.095522 -104.282887 -47.196752
37 38 39 40 41 42
-197.970080 -107.636919 -126.446916 -58.109455 -103.361898 -82.148141
43 44 45 46 47 48
-133.073015 -137.638054 -218.683608 -234.938208 -235.829191 -114.436672
49 50 51 52 53 54
-50.597597 -82.985887 -200.151222 1.280567 49.240641 131.699882
55 56 57 58 59
128.113220 211.968082 126.783891 133.080177 203.768000
> postscript(file="/var/www/html/rcomp/tmp/6tpcu1258742090.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 = 59
Frequency = 1
lag(myerror, k = 1) myerror
0 -29.219301 NA
1 -143.156225 -29.219301
2 -57.611309 -143.156225
3 -111.797698 -57.611309
4 -121.533992 -111.797698
5 -174.154439 -121.533992
6 -156.655715 -174.154439
7 -181.018954 -156.655715
8 -116.797925 -181.018954
9 7.116710 -116.797925
10 6.588924 7.116710
11 16.039076 6.588924
12 110.293749 16.039076
13 278.685619 110.293749
14 313.167184 278.685619
15 157.298001 313.167184
16 126.130083 157.298001
17 13.194445 126.130083
18 58.521328 13.194445
19 60.945005 58.521328
20 146.107288 60.945005
21 138.836843 146.107288
22 129.755155 138.836843
23 145.594349 129.755155
24 167.493228 145.594349
25 55.093412 167.493228
26 71.042262 55.093412
27 11.328585 71.042262
28 49.525166 11.328585
29 111.408253 49.525166
30 103.094182 111.408253
31 45.743921 103.094182
32 62.590354 45.743921
33 -44.095522 62.590354
34 -104.282887 -44.095522
35 -47.196752 -104.282887
36 -197.970080 -47.196752
37 -107.636919 -197.970080
38 -126.446916 -107.636919
39 -58.109455 -126.446916
40 -103.361898 -58.109455
41 -82.148141 -103.361898
42 -133.073015 -82.148141
43 -137.638054 -133.073015
44 -218.683608 -137.638054
45 -234.938208 -218.683608
46 -235.829191 -234.938208
47 -114.436672 -235.829191
48 -50.597597 -114.436672
49 -82.985887 -50.597597
50 -200.151222 -82.985887
51 1.280567 -200.151222
52 49.240641 1.280567
53 131.699882 49.240641
54 128.113220 131.699882
55 211.968082 128.113220
56 126.783891 211.968082
57 133.080177 126.783891
58 203.768000 133.080177
59 NA 203.768000
> dum1 <- dum[2:length(myerror),]
> dum1
lag(myerror, k = 1) myerror
[1,] -143.156225 -29.219301
[2,] -57.611309 -143.156225
[3,] -111.797698 -57.611309
[4,] -121.533992 -111.797698
[5,] -174.154439 -121.533992
[6,] -156.655715 -174.154439
[7,] -181.018954 -156.655715
[8,] -116.797925 -181.018954
[9,] 7.116710 -116.797925
[10,] 6.588924 7.116710
[11,] 16.039076 6.588924
[12,] 110.293749 16.039076
[13,] 278.685619 110.293749
[14,] 313.167184 278.685619
[15,] 157.298001 313.167184
[16,] 126.130083 157.298001
[17,] 13.194445 126.130083
[18,] 58.521328 13.194445
[19,] 60.945005 58.521328
[20,] 146.107288 60.945005
[21,] 138.836843 146.107288
[22,] 129.755155 138.836843
[23,] 145.594349 129.755155
[24,] 167.493228 145.594349
[25,] 55.093412 167.493228
[26,] 71.042262 55.093412
[27,] 11.328585 71.042262
[28,] 49.525166 11.328585
[29,] 111.408253 49.525166
[30,] 103.094182 111.408253
[31,] 45.743921 103.094182
[32,] 62.590354 45.743921
[33,] -44.095522 62.590354
[34,] -104.282887 -44.095522
[35,] -47.196752 -104.282887
[36,] -197.970080 -47.196752
[37,] -107.636919 -197.970080
[38,] -126.446916 -107.636919
[39,] -58.109455 -126.446916
[40,] -103.361898 -58.109455
[41,] -82.148141 -103.361898
[42,] -133.073015 -82.148141
[43,] -137.638054 -133.073015
[44,] -218.683608 -137.638054
[45,] -234.938208 -218.683608
[46,] -235.829191 -234.938208
[47,] -114.436672 -235.829191
[48,] -50.597597 -114.436672
[49,] -82.985887 -50.597597
[50,] -200.151222 -82.985887
[51,] 1.280567 -200.151222
[52,] 49.240641 1.280567
[53,] 131.699882 49.240641
[54,] 128.113220 131.699882
[55,] 211.968082 128.113220
[56,] 126.783891 211.968082
[57,] 133.080177 126.783891
[58,] 203.768000 133.080177
> z <- as.data.frame(dum1)
> z
lag(myerror, k = 1) myerror
1 -143.156225 -29.219301
2 -57.611309 -143.156225
3 -111.797698 -57.611309
4 -121.533992 -111.797698
5 -174.154439 -121.533992
6 -156.655715 -174.154439
7 -181.018954 -156.655715
8 -116.797925 -181.018954
9 7.116710 -116.797925
10 6.588924 7.116710
11 16.039076 6.588924
12 110.293749 16.039076
13 278.685619 110.293749
14 313.167184 278.685619
15 157.298001 313.167184
16 126.130083 157.298001
17 13.194445 126.130083
18 58.521328 13.194445
19 60.945005 58.521328
20 146.107288 60.945005
21 138.836843 146.107288
22 129.755155 138.836843
23 145.594349 129.755155
24 167.493228 145.594349
25 55.093412 167.493228
26 71.042262 55.093412
27 11.328585 71.042262
28 49.525166 11.328585
29 111.408253 49.525166
30 103.094182 111.408253
31 45.743921 103.094182
32 62.590354 45.743921
33 -44.095522 62.590354
34 -104.282887 -44.095522
35 -47.196752 -104.282887
36 -197.970080 -47.196752
37 -107.636919 -197.970080
38 -126.446916 -107.636919
39 -58.109455 -126.446916
40 -103.361898 -58.109455
41 -82.148141 -103.361898
42 -133.073015 -82.148141
43 -137.638054 -133.073015
44 -218.683608 -137.638054
45 -234.938208 -218.683608
46 -235.829191 -234.938208
47 -114.436672 -235.829191
48 -50.597597 -114.436672
49 -82.985887 -50.597597
50 -200.151222 -82.985887
51 1.280567 -200.151222
52 49.240641 1.280567
53 131.699882 49.240641
54 128.113220 131.699882
55 211.968082 128.113220
56 126.783891 211.968082
57 133.080177 126.783891
58 203.768000 133.080177
> 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/71kq51258742090.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/8ymst1258742090.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/94zie1258742090.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/10dpym1258742090.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/11yy0b1258742090.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/12yyjv1258742090.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/13bqci1258742090.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/14hy8t1258742090.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/15vqj31258742090.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/169kda1258742090.tab")
+ }
>
> system("convert tmp/19hva1258742090.ps tmp/19hva1258742090.png")
> system("convert tmp/2ympw1258742090.ps tmp/2ympw1258742090.png")
> system("convert tmp/3sfu31258742090.ps tmp/3sfu31258742090.png")
> system("convert tmp/4687m1258742090.ps tmp/4687m1258742090.png")
> system("convert tmp/5emfa1258742090.ps tmp/5emfa1258742090.png")
> system("convert tmp/6tpcu1258742090.ps tmp/6tpcu1258742090.png")
> system("convert tmp/71kq51258742090.ps tmp/71kq51258742090.png")
> system("convert tmp/8ymst1258742090.ps tmp/8ymst1258742090.png")
> system("convert tmp/94zie1258742090.ps tmp/94zie1258742090.png")
> system("convert tmp/10dpym1258742090.ps tmp/10dpym1258742090.png")
>
>
> proc.time()
user system elapsed
2.351 1.635 2.831