R version 2.8.0 (2008-10-20)
Copyright (C) 2008 The R Foundation for Statistical Computing
ISBN 3-900051-07-0
R is free software and comes with ABSOLUTELY NO WARRANTY.
You are welcome to redistribute it under certain conditions.
Type 'license()' or 'licence()' for distribution details.
Natural language support but running in an English locale
R is a collaborative project with many contributors.
Type 'contributors()' for more information and
'citation()' on how to cite R or R packages in publications.
Type 'demo()' for some demos, 'help()' for on-line help, or
'help.start()' for an HTML browser interface to help.
Type 'q()' to quit R.
> x <- array(list(1635.25
+ ,8169.75
+ ,7977.64
+ ,10171
+ ,-14.9
+ ,-18
+ ,1.8
+ ,2.05
+ ,1833.42
+ ,7905.84
+ ,8334.59
+ ,9721
+ ,-16.2
+ ,-11
+ ,1.5
+ ,2.05
+ ,1910.43
+ ,8145.82
+ ,8623.36
+ ,9897
+ ,-14.4
+ ,-9
+ ,1
+ ,1.81
+ ,1959.67
+ ,8895.71
+ ,9098.03
+ ,9828
+ ,-17.3
+ ,-10
+ ,1.6
+ ,1.58
+ ,1969.6
+ ,9676.31
+ ,9154.34
+ ,9924
+ ,-15.7
+ ,-13
+ ,1.5
+ ,1.57
+ ,2061.41
+ ,9884.59
+ ,9284.73
+ ,10371
+ ,-12.6
+ ,-11
+ ,1.8
+ ,1.76
+ ,2093.48
+ ,10637.44
+ ,9492.49
+ ,10846
+ ,-9.4
+ ,-5
+ ,1.8
+ ,1.76
+ ,2120.88
+ ,10717.13
+ ,9682.35
+ ,10413
+ ,-8.1
+ ,-15
+ ,1.6
+ ,1.89
+ ,2174.56
+ ,10205.29
+ ,9762.12
+ ,10709
+ ,-5.4
+ ,-6
+ ,1.9
+ ,1.9
+ ,2196.72
+ ,10295.98
+ ,10124.63
+ ,10662
+ ,-4.6
+ ,-6
+ ,1.7
+ ,1.9
+ ,2350.44
+ ,10892.76
+ ,10540.05
+ ,10570
+ ,-4.9
+ ,-3
+ ,1.6
+ ,1.92
+ ,2440.25
+ ,10631.92
+ ,10601.61
+ ,10297
+ ,-4
+ ,-1
+ ,1.3
+ ,1.76
+ ,2408.64
+ ,11441.08
+ ,10323.73
+ ,10635
+ ,-3.1
+ ,-3
+ ,1.1
+ ,1.64
+ ,2472.81
+ ,11950.95
+ ,10418.4
+ ,10872
+ ,-1.3
+ ,-4
+ ,1.9
+ ,1.57
+ ,2407.6
+ ,11037.54
+ ,10092.96
+ ,10296
+ ,0
+ ,-6
+ ,2.6
+ ,1.69
+ ,2454.62
+ ,11527.72
+ ,10364.91
+ ,10383
+ ,-0.4
+ ,0
+ ,2.3
+ ,1.76
+ ,2448.05
+ ,11383.89
+ ,10152.09
+ ,10431
+ ,3
+ ,-4
+ ,2.4
+ ,1.89
+ ,2497.84
+ ,10989.34
+ ,10032.8
+ ,10574
+ ,0.4
+ ,-2
+ ,2.2
+ ,1.78
+ ,2645.64
+ ,11079.42
+ ,10204.59
+ ,10653
+ ,1.2
+ ,-2
+ ,2
+ ,1.88
+ ,2756.76
+ ,11028.93
+ ,10001.6
+ ,10805
+ ,0.6
+ ,-6
+ ,2.9
+ ,1.86
+ ,2849.27
+ ,10973
+ ,10411.75
+ ,10872
+ ,-1.3
+ ,-7
+ ,2.6
+ ,1.88
+ ,2921.44
+ ,11068.05
+ ,10673.38
+ ,10625
+ ,-3.2
+ ,-6
+ ,2.3
+ ,1.87
+ ,2981.85
+ ,11394.84
+ ,10539.51
+ ,10407
+ ,-1.8
+ ,-6
+ ,2.3
+ ,1.86
+ ,3080.58
+ ,11545.71
+ ,10723.78
+ ,10463
+ ,-3.6
+ ,-3
+ ,2.6
+ ,1.89
+ ,3106.22
+ ,11809.38
+ ,10682.06
+ ,10556
+ ,-4.2
+ ,-2
+ ,3.1
+ ,1.9
+ ,3119.31
+ ,11395.64
+ ,10283.19
+ ,10646
+ ,-6.9
+ ,-5
+ ,2.8
+ ,1.89
+ ,3061.26
+ ,11082.38
+ ,10377.18
+ ,10702
+ ,-8
+ ,-11
+ ,2.5
+ ,1.85
+ ,3097.31
+ ,11402.75
+ ,10486.64
+ ,11353
+ ,-7.5
+ ,-11
+ ,2.9
+ ,1.78
+ ,3161.69
+ ,11716.87
+ ,10545.38
+ ,11346
+ ,-8.2
+ ,-11
+ ,3.1
+ ,1.71
+ ,3257.16
+ ,12204.98
+ ,10554.27
+ ,11451
+ ,-7.6
+ ,-10
+ ,3.1
+ ,1.69
+ ,3277.01
+ ,12986.62
+ ,10532.54
+ ,11964
+ ,-3.7
+ ,-14
+ ,3.2
+ ,1.72
+ ,3295.32
+ ,13392.79
+ ,10324.31
+ ,12574
+ ,-1.7
+ ,-8
+ ,2.5
+ ,1.77
+ ,3363.99
+ ,14368.05
+ ,10695.25
+ ,13031
+ ,-0.7
+ ,-9
+ ,2.6
+ ,1.98
+ ,3494.17
+ ,15650.83
+ ,10827.81
+ ,13812
+ ,0.2
+ ,-5
+ ,2.9
+ ,2.2
+ ,3667.03
+ ,16102.64
+ ,10872.48
+ ,14544
+ ,0.6
+ ,-1
+ ,2.6
+ ,2.25
+ ,3813.06
+ ,16187.64
+ ,10971.19
+ ,14931
+ ,2.2
+ ,-2
+ ,2.4
+ ,2.24
+ ,3917.96
+ ,16311.54
+ ,11145.65
+ ,14886
+ ,3.3
+ ,-5
+ ,1.7
+ ,2.51
+ ,3895.51
+ ,17232.97
+ ,11234.68
+ ,16005
+ ,5.3
+ ,-4
+ ,2
+ ,2.79
+ ,3801.06
+ ,16397.83
+ ,11333.88
+ ,17064
+ ,5.5
+ ,-6
+ ,2.2
+ ,3.07
+ ,3570.12
+ ,14990.31
+ ,10997.97
+ ,15168
+ ,6.3
+ ,-2
+ ,1.9
+ ,3.08
+ ,3701.61
+ ,15147.55
+ ,11036.89
+ ,16050
+ ,7.7
+ ,-2
+ ,1.6
+ ,3.05
+ ,3862.27
+ ,15786.78
+ ,11257.35
+ ,15839
+ ,6.5
+ ,-2
+ ,1.6
+ ,3.08
+ ,3970.1
+ ,15934.09
+ ,11533.59
+ ,15137
+ ,5.5
+ ,-2
+ ,1.2
+ ,3.15
+ ,4138.52
+ ,16519.44
+ ,11963.12
+ ,14954
+ ,6.9
+ ,2
+ ,1.2
+ ,3.16
+ ,4199.75
+ ,16101.07
+ ,12185.15
+ ,15648
+ ,5.7
+ ,1
+ ,1.5
+ ,3.16
+ ,4290.89
+ ,16775.08
+ ,12377.62
+ ,15305
+ ,6.9
+ ,-8
+ ,1.6
+ ,3.19
+ ,4443.91
+ ,17286.32
+ ,12512.89
+ ,15579
+ ,6.1
+ ,-1
+ ,1.7
+ ,3.44
+ ,4502.64
+ ,17741.23
+ ,12631.48
+ ,16348
+ ,4.8
+ ,1
+ ,1.8
+ ,3.55
+ ,4356.98
+ ,17128.37
+ ,12268.53
+ ,15928
+ ,3.7
+ ,-1
+ ,1.8
+ ,3.6
+ ,4591.27
+ ,17460.53
+ ,12754.8
+ ,16171
+ ,5.8
+ ,2
+ ,1.8
+ ,3.62
+ ,4696.96
+ ,17611.14
+ ,13407.75
+ ,15937
+ ,6.8
+ ,2
+ ,1.3
+ ,3.69)
+ ,dim=c(8
+ ,51)
+ ,dimnames=list(c('BEL_20'
+ ,'Nikkei'
+ ,'DJ_Indust'
+ ,'Goudprijs'
+ ,'Conjunct_Seizoenzuiver'
+ ,'Cons_vertrouw'
+ ,'Alg_consumptie_index_BE'
+ ,'Gem_rente_kasbon_1j')
+ ,1:51))
> y <- array(NA,dim=c(8,51),dimnames=list(c('BEL_20','Nikkei','DJ_Indust','Goudprijs','Conjunct_Seizoenzuiver','Cons_vertrouw','Alg_consumptie_index_BE','Gem_rente_kasbon_1j'),1:51))
> 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 = 'Do not include Seasonal 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
BEL_20 Nikkei DJ_Indust Goudprijs Conjunct_Seizoenzuiver Cons_vertrouw
1 1635.25 8169.75 7977.64 10171 -14.9 -18
2 1833.42 7905.84 8334.59 9721 -16.2 -11
3 1910.43 8145.82 8623.36 9897 -14.4 -9
4 1959.67 8895.71 9098.03 9828 -17.3 -10
5 1969.60 9676.31 9154.34 9924 -15.7 -13
6 2061.41 9884.59 9284.73 10371 -12.6 -11
7 2093.48 10637.44 9492.49 10846 -9.4 -5
8 2120.88 10717.13 9682.35 10413 -8.1 -15
9 2174.56 10205.29 9762.12 10709 -5.4 -6
10 2196.72 10295.98 10124.63 10662 -4.6 -6
11 2350.44 10892.76 10540.05 10570 -4.9 -3
12 2440.25 10631.92 10601.61 10297 -4.0 -1
13 2408.64 11441.08 10323.73 10635 -3.1 -3
14 2472.81 11950.95 10418.40 10872 -1.3 -4
15 2407.60 11037.54 10092.96 10296 0.0 -6
16 2454.62 11527.72 10364.91 10383 -0.4 0
17 2448.05 11383.89 10152.09 10431 3.0 -4
18 2497.84 10989.34 10032.80 10574 0.4 -2
19 2645.64 11079.42 10204.59 10653 1.2 -2
20 2756.76 11028.93 10001.60 10805 0.6 -6
21 2849.27 10973.00 10411.75 10872 -1.3 -7
22 2921.44 11068.05 10673.38 10625 -3.2 -6
23 2981.85 11394.84 10539.51 10407 -1.8 -6
24 3080.58 11545.71 10723.78 10463 -3.6 -3
25 3106.22 11809.38 10682.06 10556 -4.2 -2
26 3119.31 11395.64 10283.19 10646 -6.9 -5
27 3061.26 11082.38 10377.18 10702 -8.0 -11
28 3097.31 11402.75 10486.64 11353 -7.5 -11
29 3161.69 11716.87 10545.38 11346 -8.2 -11
30 3257.16 12204.98 10554.27 11451 -7.6 -10
31 3277.01 12986.62 10532.54 11964 -3.7 -14
32 3295.32 13392.79 10324.31 12574 -1.7 -8
33 3363.99 14368.05 10695.25 13031 -0.7 -9
34 3494.17 15650.83 10827.81 13812 0.2 -5
35 3667.03 16102.64 10872.48 14544 0.6 -1
36 3813.06 16187.64 10971.19 14931 2.2 -2
37 3917.96 16311.54 11145.65 14886 3.3 -5
38 3895.51 17232.97 11234.68 16005 5.3 -4
39 3801.06 16397.83 11333.88 17064 5.5 -6
40 3570.12 14990.31 10997.97 15168 6.3 -2
41 3701.61 15147.55 11036.89 16050 7.7 -2
42 3862.27 15786.78 11257.35 15839 6.5 -2
43 3970.10 15934.09 11533.59 15137 5.5 -2
44 4138.52 16519.44 11963.12 14954 6.9 2
45 4199.75 16101.07 12185.15 15648 5.7 1
46 4290.89 16775.08 12377.62 15305 6.9 -8
47 4443.91 17286.32 12512.89 15579 6.1 -1
48 4502.64 17741.23 12631.48 16348 4.8 1
49 4356.98 17128.37 12268.53 15928 3.7 -1
50 4591.27 17460.53 12754.80 16171 5.8 2
51 4696.96 17611.14 13407.75 15937 6.8 2
Alg_consumptie_index_BE Gem_rente_kasbon_1j t
1 1.8 2.05 1
2 1.5 2.05 2
3 1.0 1.81 3
4 1.6 1.58 4
5 1.5 1.57 5
6 1.8 1.76 6
7 1.8 1.76 7
8 1.6 1.89 8
9 1.9 1.90 9
10 1.7 1.90 10
11 1.6 1.92 11
12 1.3 1.76 12
13 1.1 1.64 13
14 1.9 1.57 14
15 2.6 1.69 15
16 2.3 1.76 16
17 2.4 1.89 17
18 2.2 1.78 18
19 2.0 1.88 19
20 2.9 1.86 20
21 2.6 1.88 21
22 2.3 1.87 22
23 2.3 1.86 23
24 2.6 1.89 24
25 3.1 1.90 25
26 2.8 1.89 26
27 2.5 1.85 27
28 2.9 1.78 28
29 3.1 1.71 29
30 3.1 1.69 30
31 3.2 1.72 31
32 2.5 1.77 32
33 2.6 1.98 33
34 2.9 2.20 34
35 2.6 2.25 35
36 2.4 2.24 36
37 1.7 2.51 37
38 2.0 2.79 38
39 2.2 3.07 39
40 1.9 3.08 40
41 1.6 3.05 41
42 1.6 3.08 42
43 1.2 3.15 43
44 1.2 3.16 44
45 1.5 3.16 45
46 1.6 3.19 46
47 1.7 3.44 47
48 1.8 3.55 48
49 1.8 3.60 49
50 1.8 3.62 50
51 1.3 3.69 51
> k <- length(x[1,])
> df <- as.data.frame(x)
> (mylm <- lm(df))
Call:
lm(formula = df)
Coefficients:
(Intercept) Nikkei DJ_Indust
131.46461 0.07005 0.10697
Goudprijs Conjunct_Seizoenzuiver Cons_vertrouw
-0.02675 -17.04801 4.47598
Alg_consumptie_index_BE Gem_rente_kasbon_1j t
0.15020 58.67386 43.19066
> (mysum <- summary(mylm))
Call:
lm(formula = df)
Residuals:
Min 1Q Median 3Q Max
-174.437 -46.518 -7.151 38.823 182.834
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 131.46461 438.79514 0.300 0.76596
Nikkei 0.07005 0.02418 2.898 0.00595 **
DJ_Indust 0.10697 0.03939 2.715 0.00957 **
Goudprijs -0.02675 0.02880 -0.929 0.35834
Conjunct_Seizoenzuiver -17.04801 3.71115 -4.594 3.94e-05 ***
Cons_vertrouw 4.47598 4.04205 1.107 0.27444
Alg_consumptie_index_BE 0.15020 32.40034 0.005 0.99632
Gem_rente_kasbon_1j 58.67386 58.12821 1.009 0.31857
t 43.19066 3.56291 12.122 2.66e-15 ***
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Residual standard error: 77.96 on 42 degrees of freedom
Multiple R-squared: 0.9927, Adjusted R-squared: 0.9914
F-statistic: 717.7 on 8 and 42 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.10323554 0.20647108 0.8967645
[2,] 0.07220871 0.14441742 0.9277913
[3,] 0.04493255 0.08986511 0.9550674
[4,] 0.02652611 0.05305222 0.9734739
[5,] 0.05880921 0.11761841 0.9411908
[6,] 0.04519405 0.09038809 0.9548060
[7,] 0.06299778 0.12599557 0.9370022
[8,] 0.14854519 0.29709038 0.8514548
[9,] 0.45877606 0.91755213 0.5412239
[10,] 0.36942737 0.73885473 0.6305726
[11,] 0.34698139 0.69396278 0.6530186
[12,] 0.30215190 0.60430380 0.6978481
[13,] 0.23108400 0.46216800 0.7689160
[14,] 0.16385543 0.32771085 0.8361446
[15,] 0.39586553 0.79173106 0.6041345
[16,] 0.71522966 0.56954069 0.2847703
[17,] 0.73887132 0.52225737 0.2611287
[18,] 0.67445075 0.65109851 0.3255493
[19,] 0.65791758 0.68416484 0.3420824
[20,] 0.84689729 0.30620543 0.1531027
[21,] 0.86213655 0.27572691 0.1378635
[22,] 0.78876757 0.42246486 0.2112324
[23,] 0.72910456 0.54179087 0.2708954
[24,] 0.77105591 0.45788818 0.2289441
[25,] 0.83054551 0.33890897 0.1694545
[26,] 0.87213453 0.25573094 0.1278655
[27,] 0.83613050 0.32773900 0.1638695
[28,] 0.73835928 0.52328143 0.2616407
> postscript(file="/var/www/html/freestat/rcomp/tmp/1ctmd1291647490.ps",horizontal=F,onefile=F,pagecentre=F,paper="special",width=8.3333333333333,height=5.5555555555556)
> plot(x[,1], type='l', main='Actuals and Interpolation', ylab='value of Actuals and Interpolation (dots)', xlab='time or index')
> points(x[,1]-mysum$resid)
> grid()
> dev.off()
null device
1
> postscript(file="/var/www/html/freestat/rcomp/tmp/2m2ly1291647490.ps",horizontal=F,onefile=F,pagecentre=F,paper="special",width=8.3333333333333,height=5.5555555555556)
> plot(mysum$resid, type='b', pch=19, main='Residuals', ylab='value of Residuals', xlab='time or index')
> grid()
> dev.off()
null device
1
> postscript(file="/var/www/html/freestat/rcomp/tmp/3m2ly1291647490.ps",horizontal=F,onefile=F,pagecentre=F,paper="special",width=8.3333333333333,height=5.5555555555556)
> hist(mysum$resid, main='Residual Histogram', xlab='values of Residuals')
> grid()
> dev.off()
null device
1
> postscript(file="/var/www/html/freestat/rcomp/tmp/4ft211291647490.ps",horizontal=F,onefile=F,pagecentre=F,paper="special",width=8.3333333333333,height=5.5555555555556)
> densityplot(~mysum$resid,col='black',main='Residual Density Plot', xlab='values of Residuals')
> dev.off()
null device
1
> postscript(file="/var/www/html/freestat/rcomp/tmp/5ft211291647490.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 = 51
Frequency = 1
1 2 3 4 5 6
12.992845 82.790546 109.508274 -21.154539 -71.248419 -6.506254
7 8 9 10 11 12
-52.186766 -46.127635 4.717650 -49.032379 -46.908969 19.919465
13 14 15 16 17 18
-41.433956 -20.810593 -21.851111 -116.861072 -64.272118 -60.239767
19 20 21 22 23 24
29.597273 135.556829 117.669670 69.161204 96.429676 77.267172
25 26 27 28 29 30
32.829238 44.817245 -32.539056 -43.845913 -58.989123 -32.119809
31 32 33 34 35 36
-11.554165 -21.885041 -82.993312 -94.669643 4.178444 133.226434
37 38 39 40 41 42
182.834412 86.200692 20.685920 -174.436835 -52.049294 -30.805408
43 44 45 46 47 48
-45.909083 -7.151263 18.984828 48.919207 56.140266 10.113664
49 50 51
-120.949049 22.565138 11.424481
> postscript(file="/var/www/html/freestat/rcomp/tmp/6ft211291647490.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 = 51
Frequency = 1
lag(myerror, k = 1) myerror
0 12.992845 NA
1 82.790546 12.992845
2 109.508274 82.790546
3 -21.154539 109.508274
4 -71.248419 -21.154539
5 -6.506254 -71.248419
6 -52.186766 -6.506254
7 -46.127635 -52.186766
8 4.717650 -46.127635
9 -49.032379 4.717650
10 -46.908969 -49.032379
11 19.919465 -46.908969
12 -41.433956 19.919465
13 -20.810593 -41.433956
14 -21.851111 -20.810593
15 -116.861072 -21.851111
16 -64.272118 -116.861072
17 -60.239767 -64.272118
18 29.597273 -60.239767
19 135.556829 29.597273
20 117.669670 135.556829
21 69.161204 117.669670
22 96.429676 69.161204
23 77.267172 96.429676
24 32.829238 77.267172
25 44.817245 32.829238
26 -32.539056 44.817245
27 -43.845913 -32.539056
28 -58.989123 -43.845913
29 -32.119809 -58.989123
30 -11.554165 -32.119809
31 -21.885041 -11.554165
32 -82.993312 -21.885041
33 -94.669643 -82.993312
34 4.178444 -94.669643
35 133.226434 4.178444
36 182.834412 133.226434
37 86.200692 182.834412
38 20.685920 86.200692
39 -174.436835 20.685920
40 -52.049294 -174.436835
41 -30.805408 -52.049294
42 -45.909083 -30.805408
43 -7.151263 -45.909083
44 18.984828 -7.151263
45 48.919207 18.984828
46 56.140266 48.919207
47 10.113664 56.140266
48 -120.949049 10.113664
49 22.565138 -120.949049
50 11.424481 22.565138
51 NA 11.424481
> dum1 <- dum[2:length(myerror),]
> dum1
lag(myerror, k = 1) myerror
[1,] 82.790546 12.992845
[2,] 109.508274 82.790546
[3,] -21.154539 109.508274
[4,] -71.248419 -21.154539
[5,] -6.506254 -71.248419
[6,] -52.186766 -6.506254
[7,] -46.127635 -52.186766
[8,] 4.717650 -46.127635
[9,] -49.032379 4.717650
[10,] -46.908969 -49.032379
[11,] 19.919465 -46.908969
[12,] -41.433956 19.919465
[13,] -20.810593 -41.433956
[14,] -21.851111 -20.810593
[15,] -116.861072 -21.851111
[16,] -64.272118 -116.861072
[17,] -60.239767 -64.272118
[18,] 29.597273 -60.239767
[19,] 135.556829 29.597273
[20,] 117.669670 135.556829
[21,] 69.161204 117.669670
[22,] 96.429676 69.161204
[23,] 77.267172 96.429676
[24,] 32.829238 77.267172
[25,] 44.817245 32.829238
[26,] -32.539056 44.817245
[27,] -43.845913 -32.539056
[28,] -58.989123 -43.845913
[29,] -32.119809 -58.989123
[30,] -11.554165 -32.119809
[31,] -21.885041 -11.554165
[32,] -82.993312 -21.885041
[33,] -94.669643 -82.993312
[34,] 4.178444 -94.669643
[35,] 133.226434 4.178444
[36,] 182.834412 133.226434
[37,] 86.200692 182.834412
[38,] 20.685920 86.200692
[39,] -174.436835 20.685920
[40,] -52.049294 -174.436835
[41,] -30.805408 -52.049294
[42,] -45.909083 -30.805408
[43,] -7.151263 -45.909083
[44,] 18.984828 -7.151263
[45,] 48.919207 18.984828
[46,] 56.140266 48.919207
[47,] 10.113664 56.140266
[48,] -120.949049 10.113664
[49,] 22.565138 -120.949049
[50,] 11.424481 22.565138
> z <- as.data.frame(dum1)
> z
lag(myerror, k = 1) myerror
1 82.790546 12.992845
2 109.508274 82.790546
3 -21.154539 109.508274
4 -71.248419 -21.154539
5 -6.506254 -71.248419
6 -52.186766 -6.506254
7 -46.127635 -52.186766
8 4.717650 -46.127635
9 -49.032379 4.717650
10 -46.908969 -49.032379
11 19.919465 -46.908969
12 -41.433956 19.919465
13 -20.810593 -41.433956
14 -21.851111 -20.810593
15 -116.861072 -21.851111
16 -64.272118 -116.861072
17 -60.239767 -64.272118
18 29.597273 -60.239767
19 135.556829 29.597273
20 117.669670 135.556829
21 69.161204 117.669670
22 96.429676 69.161204
23 77.267172 96.429676
24 32.829238 77.267172
25 44.817245 32.829238
26 -32.539056 44.817245
27 -43.845913 -32.539056
28 -58.989123 -43.845913
29 -32.119809 -58.989123
30 -11.554165 -32.119809
31 -21.885041 -11.554165
32 -82.993312 -21.885041
33 -94.669643 -82.993312
34 4.178444 -94.669643
35 133.226434 4.178444
36 182.834412 133.226434
37 86.200692 182.834412
38 20.685920 86.200692
39 -174.436835 20.685920
40 -52.049294 -174.436835
41 -30.805408 -52.049294
42 -45.909083 -30.805408
43 -7.151263 -45.909083
44 18.984828 -7.151263
45 48.919207 18.984828
46 56.140266 48.919207
47 10.113664 56.140266
48 -120.949049 10.113664
49 22.565138 -120.949049
50 11.424481 22.565138
> plot(z,main=paste('Residual Lag plot, lowess, and regression line'), ylab='values of Residuals', xlab='lagged values of Residuals')
> lines(lowess(z))
> abline(lm(z))
> grid()
> dev.off()
null device
1
> postscript(file="/var/www/html/freestat/rcomp/tmp/783j41291647490.ps",horizontal=F,onefile=F,pagecentre=F,paper="special",width=8.3333333333333,height=5.5555555555556)
> acf(mysum$resid, lag.max=length(mysum$resid)/2, main='Residual Autocorrelation Function')
> grid()
> dev.off()
null device
1
> postscript(file="/var/www/html/freestat/rcomp/tmp/80u1p1291647490.ps",horizontal=F,onefile=F,pagecentre=F,paper="special",width=8.3333333333333,height=5.5555555555556)
> pacf(mysum$resid, lag.max=length(mysum$resid)/2, main='Residual Partial Autocorrelation Function')
> grid()
> dev.off()
null device
1
> postscript(file="/var/www/html/freestat/rcomp/tmp/90u1p1291647490.ps",horizontal=F,onefile=F,pagecentre=F,paper="special",width=8.3333333333333,height=5.5555555555556)
> opar <- par(mfrow = c(2,2), oma = c(0, 0, 1.1, 0))
> plot(mylm, las = 1, sub='Residual Diagnostics')
> par(opar)
> dev.off()
null device
1
> if (n > n25) {
+ postscript(file="/var/www/html/freestat/rcomp/tmp/100u1p1291647490.ps",horizontal=F,onefile=F,pagecentre=F,paper="special",width=8.3333333333333,height=5.5555555555556)
+ plot(kp3:nmkm3,gqarr[,2], main='Goldfeld-Quandt test',ylab='2-sided p-value',xlab='breakpoint')
+ grid()
+ dev.off()
+ }
null device
1
>
> #Note: the /var/www/html/freestat/rcomp/createtable file can be downloaded at http://www.wessa.net/cretab
> load(file="/var/www/html/freestat/rcomp/createtable")
>
> a<-table.start()
> a<-table.row.start(a)
> a<-table.element(a, 'Multiple Linear Regression - Estimated Regression Equation', 1, TRUE)
> a<-table.row.end(a)
> myeq <- colnames(x)[1]
> myeq <- paste(myeq, '[t] = ', sep='')
> for (i in 1:k){
+ if (mysum$coefficients[i,1] > 0) myeq <- paste(myeq, '+', '')
+ myeq <- paste(myeq, mysum$coefficients[i,1], sep=' ')
+ if (rownames(mysum$coefficients)[i] != '(Intercept)') {
+ myeq <- paste(myeq, rownames(mysum$coefficients)[i], sep='')
+ if (rownames(mysum$coefficients)[i] != 't') myeq <- paste(myeq, '[t]', sep='')
+ }
+ }
> myeq <- paste(myeq, ' + e[t]')
> a<-table.row.start(a)
> a<-table.element(a, myeq)
> a<-table.row.end(a)
> a<-table.end(a)
> table.save(a,file="/var/www/html/freestat/rcomp/tmp/11e4yg1291647490.tab")
> a<-table.start()
> a<-table.row.start(a)
> a<-table.element(a,hyperlink('http://www.xycoon.com/ols1.htm','Multiple Linear Regression - Ordinary Least Squares',''), 6, TRUE)
> a<-table.row.end(a)
> a<-table.row.start(a)
> a<-table.element(a,'Variable',header=TRUE)
> a<-table.element(a,'Parameter',header=TRUE)
> a<-table.element(a,'S.D.',header=TRUE)
> a<-table.element(a,'T-STAT
H0: parameter = 0',header=TRUE)
> a<-table.element(a,'2-tail p-value',header=TRUE)
> a<-table.element(a,'1-tail p-value',header=TRUE)
> a<-table.row.end(a)
> for (i in 1:k){
+ a<-table.row.start(a)
+ a<-table.element(a,rownames(mysum$coefficients)[i],header=TRUE)
+ a<-table.element(a,mysum$coefficients[i,1])
+ a<-table.element(a, round(mysum$coefficients[i,2],6))
+ a<-table.element(a, round(mysum$coefficients[i,3],4))
+ a<-table.element(a, round(mysum$coefficients[i,4],6))
+ a<-table.element(a, round(mysum$coefficients[i,4]/2,6))
+ a<-table.row.end(a)
+ }
> a<-table.end(a)
> table.save(a,file="/var/www/html/freestat/rcomp/tmp/12pdgj1291647490.tab")
> a<-table.start()
> a<-table.row.start(a)
> a<-table.element(a, 'Multiple Linear Regression - Regression Statistics', 2, TRUE)
> a<-table.row.end(a)
> a<-table.row.start(a)
> a<-table.element(a, 'Multiple R',1,TRUE)
> a<-table.element(a, sqrt(mysum$r.squared))
> a<-table.row.end(a)
> a<-table.row.start(a)
> a<-table.element(a, 'R-squared',1,TRUE)
> a<-table.element(a, mysum$r.squared)
> a<-table.row.end(a)
> a<-table.row.start(a)
> a<-table.element(a, 'Adjusted R-squared',1,TRUE)
> a<-table.element(a, mysum$adj.r.squared)
> a<-table.row.end(a)
> a<-table.row.start(a)
> a<-table.element(a, 'F-TEST (value)',1,TRUE)
> a<-table.element(a, mysum$fstatistic[1])
> a<-table.row.end(a)
> a<-table.row.start(a)
> a<-table.element(a, 'F-TEST (DF numerator)',1,TRUE)
> a<-table.element(a, mysum$fstatistic[2])
> a<-table.row.end(a)
> a<-table.row.start(a)
> a<-table.element(a, 'F-TEST (DF denominator)',1,TRUE)
> a<-table.element(a, mysum$fstatistic[3])
> a<-table.row.end(a)
> a<-table.row.start(a)
> a<-table.element(a, 'p-value',1,TRUE)
> a<-table.element(a, 1-pf(mysum$fstatistic[1],mysum$fstatistic[2],mysum$fstatistic[3]))
> a<-table.row.end(a)
> a<-table.row.start(a)
> a<-table.element(a, 'Multiple Linear Regression - Residual Statistics', 2, TRUE)
> a<-table.row.end(a)
> a<-table.row.start(a)
> a<-table.element(a, 'Residual Standard Deviation',1,TRUE)
> a<-table.element(a, mysum$sigma)
> a<-table.row.end(a)
> a<-table.row.start(a)
> a<-table.element(a, 'Sum Squared Residuals',1,TRUE)
> a<-table.element(a, sum(myerror*myerror))
> a<-table.row.end(a)
> a<-table.end(a)
> table.save(a,file="/var/www/html/freestat/rcomp/tmp/13eevd1291647490.tab")
> a<-table.start()
> a<-table.row.start(a)
> a<-table.element(a, 'Multiple Linear Regression - Actuals, Interpolation, and Residuals', 4, TRUE)
> a<-table.row.end(a)
> a<-table.row.start(a)
> a<-table.element(a, 'Time or Index', 1, TRUE)
> a<-table.element(a, 'Actuals', 1, TRUE)
> a<-table.element(a, 'Interpolation
Forecast', 1, TRUE)
> a<-table.element(a, 'Residuals
Prediction Error', 1, TRUE)
> a<-table.row.end(a)
> for (i in 1:n) {
+ a<-table.row.start(a)
+ a<-table.element(a,i, 1, TRUE)
+ a<-table.element(a,x[i])
+ a<-table.element(a,x[i]-mysum$resid[i])
+ a<-table.element(a,mysum$resid[i])
+ a<-table.row.end(a)
+ }
> a<-table.end(a)
> table.save(a,file="/var/www/html/freestat/rcomp/tmp/14p5cf1291647490.tab")
> if (n > n25) {
+ a<-table.start()
+ a<-table.row.start(a)
+ a<-table.element(a,'Goldfeld-Quandt test for Heteroskedasticity',4,TRUE)
+ a<-table.row.end(a)
+ a<-table.row.start(a)
+ a<-table.element(a,'p-values',header=TRUE)
+ a<-table.element(a,'Alternative Hypothesis',3,header=TRUE)
+ a<-table.row.end(a)
+ a<-table.row.start(a)
+ a<-table.element(a,'breakpoint index',header=TRUE)
+ a<-table.element(a,'greater',header=TRUE)
+ a<-table.element(a,'2-sided',header=TRUE)
+ a<-table.element(a,'less',header=TRUE)
+ a<-table.row.end(a)
+ for (mypoint in kp3:nmkm3) {
+ a<-table.row.start(a)
+ a<-table.element(a,mypoint,header=TRUE)
+ a<-table.element(a,gqarr[mypoint-kp3+1,1])
+ a<-table.element(a,gqarr[mypoint-kp3+1,2])
+ a<-table.element(a,gqarr[mypoint-kp3+1,3])
+ a<-table.row.end(a)
+ }
+ a<-table.end(a)
+ table.save(a,file="/var/www/html/freestat/rcomp/tmp/15aosl1291647490.tab")
+ a<-table.start()
+ a<-table.row.start(a)
+ a<-table.element(a,'Meta Analysis of Goldfeld-Quandt test for Heteroskedasticity',4,TRUE)
+ a<-table.row.end(a)
+ a<-table.row.start(a)
+ a<-table.element(a,'Description',header=TRUE)
+ a<-table.element(a,'# significant tests',header=TRUE)
+ a<-table.element(a,'% significant tests',header=TRUE)
+ a<-table.element(a,'OK/NOK',header=TRUE)
+ a<-table.row.end(a)
+ a<-table.row.start(a)
+ a<-table.element(a,'1% type I error level',header=TRUE)
+ a<-table.element(a,numsignificant1)
+ a<-table.element(a,numsignificant1/numgqtests)
+ if (numsignificant1/numgqtests < 0.01) dum <- 'OK' else dum <- 'NOK'
+ a<-table.element(a,dum)
+ a<-table.row.end(a)
+ a<-table.row.start(a)
+ a<-table.element(a,'5% type I error level',header=TRUE)
+ a<-table.element(a,numsignificant5)
+ a<-table.element(a,numsignificant5/numgqtests)
+ if (numsignificant5/numgqtests < 0.05) dum <- 'OK' else dum <- 'NOK'
+ a<-table.element(a,dum)
+ a<-table.row.end(a)
+ a<-table.row.start(a)
+ a<-table.element(a,'10% type I error level',header=TRUE)
+ a<-table.element(a,numsignificant10)
+ a<-table.element(a,numsignificant10/numgqtests)
+ if (numsignificant10/numgqtests < 0.1) dum <- 'OK' else dum <- 'NOK'
+ a<-table.element(a,dum)
+ a<-table.row.end(a)
+ a<-table.end(a)
+ table.save(a,file="/var/www/html/freestat/rcomp/tmp/166y8c1291647490.tab")
+ }
>
> try(system("convert tmp/1ctmd1291647490.ps tmp/1ctmd1291647490.png",intern=TRUE))
character(0)
> try(system("convert tmp/2m2ly1291647490.ps tmp/2m2ly1291647490.png",intern=TRUE))
character(0)
> try(system("convert tmp/3m2ly1291647490.ps tmp/3m2ly1291647490.png",intern=TRUE))
character(0)
> try(system("convert tmp/4ft211291647490.ps tmp/4ft211291647490.png",intern=TRUE))
character(0)
> try(system("convert tmp/5ft211291647490.ps tmp/5ft211291647490.png",intern=TRUE))
character(0)
> try(system("convert tmp/6ft211291647490.ps tmp/6ft211291647490.png",intern=TRUE))
character(0)
> try(system("convert tmp/783j41291647490.ps tmp/783j41291647490.png",intern=TRUE))
character(0)
> try(system("convert tmp/80u1p1291647490.ps tmp/80u1p1291647490.png",intern=TRUE))
character(0)
> try(system("convert tmp/90u1p1291647490.ps tmp/90u1p1291647490.png",intern=TRUE))
character(0)
> try(system("convert tmp/100u1p1291647490.ps tmp/100u1p1291647490.png",intern=TRUE))
character(0)
>
>
> proc.time()
user system elapsed
3.777 2.449 4.083