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 = 'No 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
1 3016.7 2756.76 1 0 0 0 0 0 0 0 0 0 0
2 3052.4 2849.27 0 1 0 0 0 0 0 0 0 0 0
3 3099.6 2921.44 0 0 1 0 0 0 0 0 0 0 0
4 3103.3 2981.85 0 0 0 1 0 0 0 0 0 0 0
5 3119.8 3080.58 0 0 0 0 1 0 0 0 0 0 0
6 3093.7 3106.22 0 0 0 0 0 1 0 0 0 0 0
7 3164.9 3119.31 0 0 0 0 0 0 1 0 0 0 0
8 3311.5 3061.26 0 0 0 0 0 0 0 1 0 0 0
9 3410.6 3097.31 0 0 0 0 0 0 0 0 1 0 0
10 3392.6 3161.69 0 0 0 0 0 0 0 0 0 1 0
11 3338.2 3257.16 0 0 0 0 0 0 0 0 0 0 1
12 3285.1 3277.01 0 0 0 0 0 0 0 0 0 0 0
13 3294.8 3295.32 1 0 0 0 0 0 0 0 0 0 0
14 3611.2 3363.99 0 1 0 0 0 0 0 0 0 0 0
15 3611.3 3494.17 0 0 1 0 0 0 0 0 0 0 0
16 3521.0 3667.03 0 0 0 1 0 0 0 0 0 0 0
17 3519.3 3813.06 0 0 0 0 1 0 0 0 0 0 0
18 3438.3 3917.96 0 0 0 0 0 1 0 0 0 0 0
19 3534.9 3895.51 0 0 0 0 0 0 1 0 0 0 0
20 3705.8 3801.06 0 0 0 0 0 0 0 1 0 0 0
21 3807.6 3570.12 0 0 0 0 0 0 0 0 1 0 0
22 3663.0 3701.61 0 0 0 0 0 0 0 0 0 1 0
23 3604.5 3862.27 0 0 0 0 0 0 0 0 0 0 1
24 3563.8 3970.10 0 0 0 0 0 0 0 0 0 0 0
25 3511.4 4138.52 1 0 0 0 0 0 0 0 0 0 0
26 3546.5 4199.75 0 1 0 0 0 0 0 0 0 0 0
27 3525.4 4290.89 0 0 1 0 0 0 0 0 0 0 0
28 3529.9 4443.91 0 0 0 1 0 0 0 0 0 0 0
29 3591.6 4502.64 0 0 0 0 1 0 0 0 0 0 0
30 3668.3 4356.98 0 0 0 0 0 1 0 0 0 0 0
31 3728.8 4591.27 0 0 0 0 0 0 1 0 0 0 0
32 3853.6 4696.96 0 0 0 0 0 0 0 1 0 0 0
33 3897.7 4621.40 0 0 0 0 0 0 0 0 1 0 0
34 3640.7 4562.84 0 0 0 0 0 0 0 0 0 1 0
35 3495.5 4202.52 0 0 0 0 0 0 0 0 0 0 1
36 3495.1 4296.49 0 0 0 0 0 0 0 0 0 0 0
37 3268.0 4435.23 1 0 0 0 0 0 0 0 0 0 0
38 3479.1 4105.18 0 1 0 0 0 0 0 0 0 0 0
39 3417.8 4116.68 0 0 1 0 0 0 0 0 0 0 0
40 3521.3 3844.49 0 0 0 1 0 0 0 0 0 0 0
41 3487.1 3720.98 0 0 0 0 1 0 0 0 0 0 0
42 3529.9 3674.40 0 0 0 0 0 1 0 0 0 0 0
43 3544.3 3857.62 0 0 0 0 0 0 1 0 0 0 0
44 3710.8 3801.06 0 0 0 0 0 0 0 1 0 0 0
45 3641.9 3504.37 0 0 0 0 0 0 0 0 1 0 0
46 3447.1 3032.60 0 0 0 0 0 0 0 0 0 1 0
47 3386.8 3047.03 0 0 0 0 0 0 0 0 0 0 1
48 3438.5 2962.34 0 0 0 0 0 0 0 0 0 0 0
49 3364.3 2197.82 1 0 0 0 0 0 0 0 0 0 0
50 3462.7 2014.45 0 1 0 0 0 0 0 0 0 0 0
51 3291.9 1862.83 0 0 1 0 0 0 0 0 0 0 0
52 3550.0 1905.41 0 0 0 1 0 0 0 0 0 0 0
53 3611.0 1810.99 0 0 0 0 1 0 0 0 0 0 0
54 3708.6 1670.07 0 0 0 0 0 1 0 0 0 0 0
55 3771.1 1864.44 0 0 0 0 0 0 1 0 0 0 0
56 4042.7 2052.02 0 0 0 0 0 0 0 1 0 0 0
57 3988.4 2029.60 0 0 0 0 0 0 0 0 1 0 0
58 3851.2 2070.83 0 0 0 0 0 0 0 0 0 1 0
59 3876.7 2293.41 0 0 0 0 0 0 0 0 0 0 1
> k <- length(x[1,])
> df <- as.data.frame(x)
> (mylm <- lm(df))
Call:
lm(formula = df)
Coefficients:
(Intercept) `Bel20\r` M1 M2 M3 M4
3401.75545 0.01210 -151.41855 -11.37448 -52.92555 2.59538
M5 M6 M7 M8 M9 M10
23.04838 45.53860 105.12086 280.99712 306.78350 157.17294
M11
98.27160
> (mysum <- summary(mylm))
Call:
lm(formula = df)
Residuals:
Min 1Q Median 3Q Max
-408.28 -94.25 39.06 118.31 348.93
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 3401.75545 161.66159 21.042 <2e-16 ***
`Bel20\r` 0.01210 0.03337 0.362 0.7186
M1 -151.41855 144.06916 -1.051 0.2987
M2 -11.37448 144.19996 -0.079 0.9375
M3 -52.92555 144.12778 -0.367 0.7151
M4 2.59538 144.06151 0.018 0.9857
M5 23.04838 144.02852 0.160 0.8736
M6 45.53860 144.11031 0.316 0.7534
M7 105.12086 143.90426 0.730 0.4688
M8 280.99712 143.88439 1.953 0.0569 .
M9 306.78350 144.06951 2.129 0.0386 *
M10 157.17294 144.20147 1.090 0.2814
M11 98.27160 144.13842 0.682 0.4988
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 214.4 on 46 degrees of freedom
Multiple R-squared: 0.3143, Adjusted R-squared: 0.1354
F-statistic: 1.757 on 12 and 46 DF, p-value: 0.08518
> 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.39491419 0.7898284 0.6050858
[2,] 0.31569457 0.6313891 0.6843054
[3,] 0.35434885 0.7086977 0.6456511
[4,] 0.28040248 0.5608050 0.7195975
[5,] 0.19844469 0.3968894 0.8015553
[6,] 0.15362627 0.3072525 0.8463737
[7,] 0.09643829 0.1928766 0.9035617
[8,] 0.06478879 0.1295776 0.9352112
[9,] 0.05082853 0.1016571 0.9491715
[10,] 0.12599889 0.2519978 0.8740011
[11,] 0.26950391 0.5390078 0.7304961
[12,] 0.41977619 0.8395524 0.5802238
[13,] 0.41680491 0.8336098 0.5831951
[14,] 0.36112477 0.7222495 0.6388752
[15,] 0.30717120 0.6143424 0.6928288
[16,] 0.25581916 0.5116383 0.7441808
[17,] 0.20250058 0.4050012 0.7974994
[18,] 0.20825290 0.4165058 0.7917471
[19,] 0.23615374 0.4723075 0.7638463
[20,] 0.19470040 0.3894008 0.8052996
[21,] 0.16493197 0.3298639 0.8350680
[22,] 0.18817875 0.3763575 0.8118213
[23,] 0.16392003 0.3278401 0.8360800
[24,] 0.28481209 0.5696242 0.7151879
[25,] 0.30338574 0.6067715 0.6966143
[26,] 0.27791948 0.5558390 0.7220805
[27,] 0.27983038 0.5596608 0.7201696
[28,] 0.32446316 0.6489263 0.6755368
> postscript(file="/var/www/html/rcomp/tmp/1osat1258741484.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/270dv1258741484.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/31jsl1258741484.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/4gpas1258741484.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/5ugfm1258741484.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
-266.9853932 -372.4485547 -284.5705314 -337.1222392 -342.2695713 -391.1699566
7 8 9 10 11 12
-379.7105705 -408.2846028 -335.4070795 -204.5753238 -201.2288789 -156.2974045
13 14 15 16 17 18
4.5996521 180.1248828 220.2011597 72.2891453 48.3696256 -56.3895670
19 20 21 22 23 24
-19.1002539 -22.9339558 55.8733428 59.2932696 57.7511116 114.0182929
25 26 27 28 29 30
210.9994703 105.3147026 124.6632461 71.7912360 112.3277834 168.2996126
31 32 33 34 35 36
166.3831445 114.0283511 133.2560194 26.5749791 -55.3648891 41.3699564
37 38 39 40 41 42
-35.9898275 39.0587150 19.1706626 70.4424135 17.2835164 38.1567758
43 44 45 46 47 48
-9.2418990 -17.9339558 -109.0312801 -148.5137234 -150.0869385 0.9091552
49 50 51 52 53 54
87.3760983 47.9502543 -79.4645371 122.5994444 164.2886458 241.1031352
55 56 57 58 59
241.6695790 335.1241631 255.3089974 267.2207985 348.9295949
> postscript(file="/var/www/html/rcomp/tmp/67sza1258741484.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 -266.9853932 NA
1 -372.4485547 -266.9853932
2 -284.5705314 -372.4485547
3 -337.1222392 -284.5705314
4 -342.2695713 -337.1222392
5 -391.1699566 -342.2695713
6 -379.7105705 -391.1699566
7 -408.2846028 -379.7105705
8 -335.4070795 -408.2846028
9 -204.5753238 -335.4070795
10 -201.2288789 -204.5753238
11 -156.2974045 -201.2288789
12 4.5996521 -156.2974045
13 180.1248828 4.5996521
14 220.2011597 180.1248828
15 72.2891453 220.2011597
16 48.3696256 72.2891453
17 -56.3895670 48.3696256
18 -19.1002539 -56.3895670
19 -22.9339558 -19.1002539
20 55.8733428 -22.9339558
21 59.2932696 55.8733428
22 57.7511116 59.2932696
23 114.0182929 57.7511116
24 210.9994703 114.0182929
25 105.3147026 210.9994703
26 124.6632461 105.3147026
27 71.7912360 124.6632461
28 112.3277834 71.7912360
29 168.2996126 112.3277834
30 166.3831445 168.2996126
31 114.0283511 166.3831445
32 133.2560194 114.0283511
33 26.5749791 133.2560194
34 -55.3648891 26.5749791
35 41.3699564 -55.3648891
36 -35.9898275 41.3699564
37 39.0587150 -35.9898275
38 19.1706626 39.0587150
39 70.4424135 19.1706626
40 17.2835164 70.4424135
41 38.1567758 17.2835164
42 -9.2418990 38.1567758
43 -17.9339558 -9.2418990
44 -109.0312801 -17.9339558
45 -148.5137234 -109.0312801
46 -150.0869385 -148.5137234
47 0.9091552 -150.0869385
48 87.3760983 0.9091552
49 47.9502543 87.3760983
50 -79.4645371 47.9502543
51 122.5994444 -79.4645371
52 164.2886458 122.5994444
53 241.1031352 164.2886458
54 241.6695790 241.1031352
55 335.1241631 241.6695790
56 255.3089974 335.1241631
57 267.2207985 255.3089974
58 348.9295949 267.2207985
59 NA 348.9295949
> dum1 <- dum[2:length(myerror),]
> dum1
lag(myerror, k = 1) myerror
[1,] -372.4485547 -266.9853932
[2,] -284.5705314 -372.4485547
[3,] -337.1222392 -284.5705314
[4,] -342.2695713 -337.1222392
[5,] -391.1699566 -342.2695713
[6,] -379.7105705 -391.1699566
[7,] -408.2846028 -379.7105705
[8,] -335.4070795 -408.2846028
[9,] -204.5753238 -335.4070795
[10,] -201.2288789 -204.5753238
[11,] -156.2974045 -201.2288789
[12,] 4.5996521 -156.2974045
[13,] 180.1248828 4.5996521
[14,] 220.2011597 180.1248828
[15,] 72.2891453 220.2011597
[16,] 48.3696256 72.2891453
[17,] -56.3895670 48.3696256
[18,] -19.1002539 -56.3895670
[19,] -22.9339558 -19.1002539
[20,] 55.8733428 -22.9339558
[21,] 59.2932696 55.8733428
[22,] 57.7511116 59.2932696
[23,] 114.0182929 57.7511116
[24,] 210.9994703 114.0182929
[25,] 105.3147026 210.9994703
[26,] 124.6632461 105.3147026
[27,] 71.7912360 124.6632461
[28,] 112.3277834 71.7912360
[29,] 168.2996126 112.3277834
[30,] 166.3831445 168.2996126
[31,] 114.0283511 166.3831445
[32,] 133.2560194 114.0283511
[33,] 26.5749791 133.2560194
[34,] -55.3648891 26.5749791
[35,] 41.3699564 -55.3648891
[36,] -35.9898275 41.3699564
[37,] 39.0587150 -35.9898275
[38,] 19.1706626 39.0587150
[39,] 70.4424135 19.1706626
[40,] 17.2835164 70.4424135
[41,] 38.1567758 17.2835164
[42,] -9.2418990 38.1567758
[43,] -17.9339558 -9.2418990
[44,] -109.0312801 -17.9339558
[45,] -148.5137234 -109.0312801
[46,] -150.0869385 -148.5137234
[47,] 0.9091552 -150.0869385
[48,] 87.3760983 0.9091552
[49,] 47.9502543 87.3760983
[50,] -79.4645371 47.9502543
[51,] 122.5994444 -79.4645371
[52,] 164.2886458 122.5994444
[53,] 241.1031352 164.2886458
[54,] 241.6695790 241.1031352
[55,] 335.1241631 241.6695790
[56,] 255.3089974 335.1241631
[57,] 267.2207985 255.3089974
[58,] 348.9295949 267.2207985
> z <- as.data.frame(dum1)
> z
lag(myerror, k = 1) myerror
1 -372.4485547 -266.9853932
2 -284.5705314 -372.4485547
3 -337.1222392 -284.5705314
4 -342.2695713 -337.1222392
5 -391.1699566 -342.2695713
6 -379.7105705 -391.1699566
7 -408.2846028 -379.7105705
8 -335.4070795 -408.2846028
9 -204.5753238 -335.4070795
10 -201.2288789 -204.5753238
11 -156.2974045 -201.2288789
12 4.5996521 -156.2974045
13 180.1248828 4.5996521
14 220.2011597 180.1248828
15 72.2891453 220.2011597
16 48.3696256 72.2891453
17 -56.3895670 48.3696256
18 -19.1002539 -56.3895670
19 -22.9339558 -19.1002539
20 55.8733428 -22.9339558
21 59.2932696 55.8733428
22 57.7511116 59.2932696
23 114.0182929 57.7511116
24 210.9994703 114.0182929
25 105.3147026 210.9994703
26 124.6632461 105.3147026
27 71.7912360 124.6632461
28 112.3277834 71.7912360
29 168.2996126 112.3277834
30 166.3831445 168.2996126
31 114.0283511 166.3831445
32 133.2560194 114.0283511
33 26.5749791 133.2560194
34 -55.3648891 26.5749791
35 41.3699564 -55.3648891
36 -35.9898275 41.3699564
37 39.0587150 -35.9898275
38 19.1706626 39.0587150
39 70.4424135 19.1706626
40 17.2835164 70.4424135
41 38.1567758 17.2835164
42 -9.2418990 38.1567758
43 -17.9339558 -9.2418990
44 -109.0312801 -17.9339558
45 -148.5137234 -109.0312801
46 -150.0869385 -148.5137234
47 0.9091552 -150.0869385
48 87.3760983 0.9091552
49 47.9502543 87.3760983
50 -79.4645371 47.9502543
51 122.5994444 -79.4645371
52 164.2886458 122.5994444
53 241.1031352 164.2886458
54 241.6695790 241.1031352
55 335.1241631 241.6695790
56 255.3089974 335.1241631
57 267.2207985 255.3089974
58 348.9295949 267.2207985
> 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/71o3u1258741484.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/8odkh1258741484.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/9va7n1258741484.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/10ox3y1258741484.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/11kx5o1258741484.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/124up31258741484.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/13oohv1258741484.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/14xgsx1258741484.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/150ok01258741484.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/1635e81258741484.tab")
+ }
>
> system("convert tmp/1osat1258741484.ps tmp/1osat1258741484.png")
> system("convert tmp/270dv1258741484.ps tmp/270dv1258741484.png")
> system("convert tmp/31jsl1258741484.ps tmp/31jsl1258741484.png")
> system("convert tmp/4gpas1258741484.ps tmp/4gpas1258741484.png")
> system("convert tmp/5ugfm1258741484.ps tmp/5ugfm1258741484.png")
> system("convert tmp/67sza1258741484.ps tmp/67sza1258741484.png")
> system("convert tmp/71o3u1258741484.ps tmp/71o3u1258741484.png")
> system("convert tmp/8odkh1258741484.ps tmp/8odkh1258741484.png")
> system("convert tmp/9va7n1258741484.ps tmp/9va7n1258741484.png")
> system("convert tmp/10ox3y1258741484.ps tmp/10ox3y1258741484.png")
>
>
> proc.time()
user system elapsed
2.368 1.552 2.784