R version 2.7.0 (2008-04-22)
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.
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(103.3,0,101.2,0,107.7,0,110.4,0,101.9,0,115.9,0,89.9,0,88.6,0,117.2,0,123.9,0,100,0,103.6,0,94.1,0,98.7,0,119.5,0,112.7,0,104.4,0,124.7,0,89.1,0,97,0,121.6,0,118.8,0,114,0,111.5,0,97.2,0,102.5,0,113.4,0,109.8,0,104.9,0,126.1,0,80,0,96.8,0,117.2,1,112.3,1,117.3,1,111.1,1,102.2,1,104.3,1,122.9,1,107.6,1,121.3,1,131.5,1,89,1,104.4,1,128.9,1,135.9,1,133.3,1,121.3,1,120.5,1,120.4,1,137.9,1,126.1,1,133.2,1,151.1,1,105,1,119,1,140.4,1,156.6,1,137.1,1,122.7,1),dim=c(2,60),dimnames=list(c('metaal','conjunctuur'),1:60))
> y <- array(NA,dim=c(2,60),dimnames=list(c('metaal','conjunctuur'),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 = '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
metaal conjunctuur M1 M2 M3 M4 M5 M6 M7 M8 M9 M10 M11 t
1 103.3 0 1 0 0 0 0 0 0 0 0 0 0 1
2 101.2 0 0 1 0 0 0 0 0 0 0 0 0 2
3 107.7 0 0 0 1 0 0 0 0 0 0 0 0 3
4 110.4 0 0 0 0 1 0 0 0 0 0 0 0 4
5 101.9 0 0 0 0 0 1 0 0 0 0 0 0 5
6 115.9 0 0 0 0 0 0 1 0 0 0 0 0 6
7 89.9 0 0 0 0 0 0 0 1 0 0 0 0 7
8 88.6 0 0 0 0 0 0 0 0 1 0 0 0 8
9 117.2 0 0 0 0 0 0 0 0 0 1 0 0 9
10 123.9 0 0 0 0 0 0 0 0 0 0 1 0 10
11 100.0 0 0 0 0 0 0 0 0 0 0 0 1 11
12 103.6 0 0 0 0 0 0 0 0 0 0 0 0 12
13 94.1 0 1 0 0 0 0 0 0 0 0 0 0 13
14 98.7 0 0 1 0 0 0 0 0 0 0 0 0 14
15 119.5 0 0 0 1 0 0 0 0 0 0 0 0 15
16 112.7 0 0 0 0 1 0 0 0 0 0 0 0 16
17 104.4 0 0 0 0 0 1 0 0 0 0 0 0 17
18 124.7 0 0 0 0 0 0 1 0 0 0 0 0 18
19 89.1 0 0 0 0 0 0 0 1 0 0 0 0 19
20 97.0 0 0 0 0 0 0 0 0 1 0 0 0 20
21 121.6 0 0 0 0 0 0 0 0 0 1 0 0 21
22 118.8 0 0 0 0 0 0 0 0 0 0 1 0 22
23 114.0 0 0 0 0 0 0 0 0 0 0 0 1 23
24 111.5 0 0 0 0 0 0 0 0 0 0 0 0 24
25 97.2 0 1 0 0 0 0 0 0 0 0 0 0 25
26 102.5 0 0 1 0 0 0 0 0 0 0 0 0 26
27 113.4 0 0 0 1 0 0 0 0 0 0 0 0 27
28 109.8 0 0 0 0 1 0 0 0 0 0 0 0 28
29 104.9 0 0 0 0 0 1 0 0 0 0 0 0 29
30 126.1 0 0 0 0 0 0 1 0 0 0 0 0 30
31 80.0 0 0 0 0 0 0 0 1 0 0 0 0 31
32 96.8 0 0 0 0 0 0 0 0 1 0 0 0 32
33 117.2 1 0 0 0 0 0 0 0 0 1 0 0 33
34 112.3 1 0 0 0 0 0 0 0 0 0 1 0 34
35 117.3 1 0 0 0 0 0 0 0 0 0 0 1 35
36 111.1 1 0 0 0 0 0 0 0 0 0 0 0 36
37 102.2 1 1 0 0 0 0 0 0 0 0 0 0 37
38 104.3 1 0 1 0 0 0 0 0 0 0 0 0 38
39 122.9 1 0 0 1 0 0 0 0 0 0 0 0 39
40 107.6 1 0 0 0 1 0 0 0 0 0 0 0 40
41 121.3 1 0 0 0 0 1 0 0 0 0 0 0 41
42 131.5 1 0 0 0 0 0 1 0 0 0 0 0 42
43 89.0 1 0 0 0 0 0 0 1 0 0 0 0 43
44 104.4 1 0 0 0 0 0 0 0 1 0 0 0 44
45 128.9 1 0 0 0 0 0 0 0 0 1 0 0 45
46 135.9 1 0 0 0 0 0 0 0 0 0 1 0 46
47 133.3 1 0 0 0 0 0 0 0 0 0 0 1 47
48 121.3 1 0 0 0 0 0 0 0 0 0 0 0 48
49 120.5 1 1 0 0 0 0 0 0 0 0 0 0 49
50 120.4 1 0 1 0 0 0 0 0 0 0 0 0 50
51 137.9 1 0 0 1 0 0 0 0 0 0 0 0 51
52 126.1 1 0 0 0 1 0 0 0 0 0 0 0 52
53 133.2 1 0 0 0 0 1 0 0 0 0 0 0 53
54 151.1 1 0 0 0 0 0 1 0 0 0 0 0 54
55 105.0 1 0 0 0 0 0 0 1 0 0 0 0 55
56 119.0 1 0 0 0 0 0 0 0 1 0 0 0 56
57 140.4 1 0 0 0 0 0 0 0 0 1 0 0 57
58 156.6 1 0 0 0 0 0 0 0 0 0 1 0 58
59 137.1 1 0 0 0 0 0 0 0 0 0 0 1 59
60 122.7 1 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) conjunctuur M1 M2 M3 M4
96.0083 -0.9972 -5.0869 -3.6444 10.6981 3.2206
M5 M6 M7 M8 M9 M10
2.5231 18.7256 -21.0519 -11.0094 12.5725 16.4950
M11 t
6.8175 0.5175
> (mysum <- summary(mylm))
Call:
lm(formula = df)
Residuals:
Min 1Q Median 3Q Max
-16.801 -3.932 0.025 4.418 15.079
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 96.0083 3.8764 24.767 < 2e-16 ***
conjunctuur -0.9972 3.7301 -0.267 0.790397
M1 -5.0869 4.5237 -1.124 0.266635
M2 -3.6444 4.5122 -0.808 0.423427
M3 10.6981 4.5032 2.376 0.021740 *
M4 3.2206 4.4968 0.716 0.477491
M5 2.5231 4.4929 0.562 0.577137
M6 18.7256 4.4916 4.169 0.000134 ***
M7 -21.0519 4.4929 -4.686 2.50e-05 ***
M8 -11.0094 4.4968 -2.448 0.018226 *
M9 12.5725 4.4877 2.802 0.007416 **
M10 16.4950 4.4813 3.681 0.000609 ***
M11 6.8175 4.4774 1.523 0.134690
t 0.5175 0.1077 4.806 1.68e-05 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 7.077 on 46 degrees of freedom
Multiple R-squared: 0.8407, Adjusted R-squared: 0.7956
F-statistic: 18.67 on 13 and 46 DF, p-value: 3.983e-14
> 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.6289741 0.7420517 0.3710259
[2,] 0.5676008 0.8647983 0.4323992
[3,] 0.5507881 0.8984238 0.4492119
[4,] 0.5394916 0.9210168 0.4605084
[5,] 0.4942793 0.9885585 0.5057207
[6,] 0.4786482 0.9572964 0.5213518
[7,] 0.5708133 0.8583734 0.4291867
[8,] 0.6635132 0.6729736 0.3364868
[9,] 0.6234482 0.7531036 0.3765518
[10,] 0.5510267 0.8979465 0.4489733
[11,] 0.4708770 0.9417540 0.5291230
[12,] 0.5006126 0.9987747 0.4993874
[13,] 0.4414024 0.8828048 0.5585976
[14,] 0.3584562 0.7169124 0.6415438
[15,] 0.4494759 0.8989518 0.5505241
[16,] 0.3560688 0.7121376 0.6439312
[17,] 0.2768268 0.5536535 0.7231732
[18,] 0.5107318 0.9785363 0.4892682
[19,] 0.5755145 0.8489710 0.4244855
[20,] 0.6432298 0.7135404 0.3567702
[21,] 0.5740160 0.8519680 0.4259840
[22,] 0.4715241 0.9430482 0.5284759
[23,] 0.3940359 0.7880717 0.6059641
[24,] 0.3580764 0.7161528 0.6419236
[25,] 0.3539006 0.7078012 0.6460994
[26,] 0.3271412 0.6542825 0.6728588
[27,] 0.2323456 0.4646913 0.7676544
> postscript(file="/var/www/html/rcomp/tmp/1w4hg1229032639.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/2gncz1229032639.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/3bp1h1229032639.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/4vllc1229032639.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/5ijr41229032639.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
11.8611111 7.8011111 -0.5588889 9.1011111 0.7811111 -1.9388889
7 8 9 10 11 12
11.3211111 -0.5388889 3.9616667 6.2216667 -8.5183333 1.3816667
13 14 15 16 17 18
-3.5488889 -0.9088889 5.0311111 5.1911111 -2.9288889 0.6511111
19 20 21 22 23 24
4.3111111 1.6511111 2.1516667 -5.0883333 -0.7283333 3.0716667
25 26 27 28 29 30
-6.6588889 -3.3188889 -7.2788889 -3.9188889 -8.6388889 -4.1588889
31 32 33 34 35 36
-10.9988889 -4.7588889 -7.4611111 -16.8011111 -2.6411111 -2.5411111
37 38 39 40 41 42
-6.8716667 -6.7316667 -2.9916667 -11.3316667 2.5483333 -3.9716667
43 44 45 46 47 48
-7.2116667 -2.3716667 -1.9711111 0.5888889 7.1488889 1.4488889
49 50 51 52 53 54
5.2183333 3.1583333 5.7983333 0.9583333 8.2383333 9.4183333
55 56 57 58 59 60
2.5783333 6.0183333 3.3188889 15.0788889 4.7388889 -3.3611111
> postscript(file="/var/www/html/rcomp/tmp/6nk4e1229032639.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 11.8611111 NA
1 7.8011111 11.8611111
2 -0.5588889 7.8011111
3 9.1011111 -0.5588889
4 0.7811111 9.1011111
5 -1.9388889 0.7811111
6 11.3211111 -1.9388889
7 -0.5388889 11.3211111
8 3.9616667 -0.5388889
9 6.2216667 3.9616667
10 -8.5183333 6.2216667
11 1.3816667 -8.5183333
12 -3.5488889 1.3816667
13 -0.9088889 -3.5488889
14 5.0311111 -0.9088889
15 5.1911111 5.0311111
16 -2.9288889 5.1911111
17 0.6511111 -2.9288889
18 4.3111111 0.6511111
19 1.6511111 4.3111111
20 2.1516667 1.6511111
21 -5.0883333 2.1516667
22 -0.7283333 -5.0883333
23 3.0716667 -0.7283333
24 -6.6588889 3.0716667
25 -3.3188889 -6.6588889
26 -7.2788889 -3.3188889
27 -3.9188889 -7.2788889
28 -8.6388889 -3.9188889
29 -4.1588889 -8.6388889
30 -10.9988889 -4.1588889
31 -4.7588889 -10.9988889
32 -7.4611111 -4.7588889
33 -16.8011111 -7.4611111
34 -2.6411111 -16.8011111
35 -2.5411111 -2.6411111
36 -6.8716667 -2.5411111
37 -6.7316667 -6.8716667
38 -2.9916667 -6.7316667
39 -11.3316667 -2.9916667
40 2.5483333 -11.3316667
41 -3.9716667 2.5483333
42 -7.2116667 -3.9716667
43 -2.3716667 -7.2116667
44 -1.9711111 -2.3716667
45 0.5888889 -1.9711111
46 7.1488889 0.5888889
47 1.4488889 7.1488889
48 5.2183333 1.4488889
49 3.1583333 5.2183333
50 5.7983333 3.1583333
51 0.9583333 5.7983333
52 8.2383333 0.9583333
53 9.4183333 8.2383333
54 2.5783333 9.4183333
55 6.0183333 2.5783333
56 3.3188889 6.0183333
57 15.0788889 3.3188889
58 4.7388889 15.0788889
59 -3.3611111 4.7388889
60 NA -3.3611111
> dum1 <- dum[2:length(myerror),]
> dum1
lag(myerror, k = 1) myerror
[1,] 7.8011111 11.8611111
[2,] -0.5588889 7.8011111
[3,] 9.1011111 -0.5588889
[4,] 0.7811111 9.1011111
[5,] -1.9388889 0.7811111
[6,] 11.3211111 -1.9388889
[7,] -0.5388889 11.3211111
[8,] 3.9616667 -0.5388889
[9,] 6.2216667 3.9616667
[10,] -8.5183333 6.2216667
[11,] 1.3816667 -8.5183333
[12,] -3.5488889 1.3816667
[13,] -0.9088889 -3.5488889
[14,] 5.0311111 -0.9088889
[15,] 5.1911111 5.0311111
[16,] -2.9288889 5.1911111
[17,] 0.6511111 -2.9288889
[18,] 4.3111111 0.6511111
[19,] 1.6511111 4.3111111
[20,] 2.1516667 1.6511111
[21,] -5.0883333 2.1516667
[22,] -0.7283333 -5.0883333
[23,] 3.0716667 -0.7283333
[24,] -6.6588889 3.0716667
[25,] -3.3188889 -6.6588889
[26,] -7.2788889 -3.3188889
[27,] -3.9188889 -7.2788889
[28,] -8.6388889 -3.9188889
[29,] -4.1588889 -8.6388889
[30,] -10.9988889 -4.1588889
[31,] -4.7588889 -10.9988889
[32,] -7.4611111 -4.7588889
[33,] -16.8011111 -7.4611111
[34,] -2.6411111 -16.8011111
[35,] -2.5411111 -2.6411111
[36,] -6.8716667 -2.5411111
[37,] -6.7316667 -6.8716667
[38,] -2.9916667 -6.7316667
[39,] -11.3316667 -2.9916667
[40,] 2.5483333 -11.3316667
[41,] -3.9716667 2.5483333
[42,] -7.2116667 -3.9716667
[43,] -2.3716667 -7.2116667
[44,] -1.9711111 -2.3716667
[45,] 0.5888889 -1.9711111
[46,] 7.1488889 0.5888889
[47,] 1.4488889 7.1488889
[48,] 5.2183333 1.4488889
[49,] 3.1583333 5.2183333
[50,] 5.7983333 3.1583333
[51,] 0.9583333 5.7983333
[52,] 8.2383333 0.9583333
[53,] 9.4183333 8.2383333
[54,] 2.5783333 9.4183333
[55,] 6.0183333 2.5783333
[56,] 3.3188889 6.0183333
[57,] 15.0788889 3.3188889
[58,] 4.7388889 15.0788889
[59,] -3.3611111 4.7388889
> z <- as.data.frame(dum1)
> z
lag(myerror, k = 1) myerror
1 7.8011111 11.8611111
2 -0.5588889 7.8011111
3 9.1011111 -0.5588889
4 0.7811111 9.1011111
5 -1.9388889 0.7811111
6 11.3211111 -1.9388889
7 -0.5388889 11.3211111
8 3.9616667 -0.5388889
9 6.2216667 3.9616667
10 -8.5183333 6.2216667
11 1.3816667 -8.5183333
12 -3.5488889 1.3816667
13 -0.9088889 -3.5488889
14 5.0311111 -0.9088889
15 5.1911111 5.0311111
16 -2.9288889 5.1911111
17 0.6511111 -2.9288889
18 4.3111111 0.6511111
19 1.6511111 4.3111111
20 2.1516667 1.6511111
21 -5.0883333 2.1516667
22 -0.7283333 -5.0883333
23 3.0716667 -0.7283333
24 -6.6588889 3.0716667
25 -3.3188889 -6.6588889
26 -7.2788889 -3.3188889
27 -3.9188889 -7.2788889
28 -8.6388889 -3.9188889
29 -4.1588889 -8.6388889
30 -10.9988889 -4.1588889
31 -4.7588889 -10.9988889
32 -7.4611111 -4.7588889
33 -16.8011111 -7.4611111
34 -2.6411111 -16.8011111
35 -2.5411111 -2.6411111
36 -6.8716667 -2.5411111
37 -6.7316667 -6.8716667
38 -2.9916667 -6.7316667
39 -11.3316667 -2.9916667
40 2.5483333 -11.3316667
41 -3.9716667 2.5483333
42 -7.2116667 -3.9716667
43 -2.3716667 -7.2116667
44 -1.9711111 -2.3716667
45 0.5888889 -1.9711111
46 7.1488889 0.5888889
47 1.4488889 7.1488889
48 5.2183333 1.4488889
49 3.1583333 5.2183333
50 5.7983333 3.1583333
51 0.9583333 5.7983333
52 8.2383333 0.9583333
53 9.4183333 8.2383333
54 2.5783333 9.4183333
55 6.0183333 2.5783333
56 3.3188889 6.0183333
57 15.0788889 3.3188889
58 4.7388889 15.0788889
59 -3.3611111 4.7388889
> 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/7bxq61229032639.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/84lfi1229032639.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/942061229032639.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/10reu01229032639.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/112ai91229032639.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/12k7ze1229032639.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/13qj7u1229032640.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/144c2k1229032640.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/151qdn1229032640.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/16e7011229032640.tab")
+ }
>
> system("convert tmp/1w4hg1229032639.ps tmp/1w4hg1229032639.png")
> system("convert tmp/2gncz1229032639.ps tmp/2gncz1229032639.png")
> system("convert tmp/3bp1h1229032639.ps tmp/3bp1h1229032639.png")
> system("convert tmp/4vllc1229032639.ps tmp/4vllc1229032639.png")
> system("convert tmp/5ijr41229032639.ps tmp/5ijr41229032639.png")
> system("convert tmp/6nk4e1229032639.ps tmp/6nk4e1229032639.png")
> system("convert tmp/7bxq61229032639.ps tmp/7bxq61229032639.png")
> system("convert tmp/84lfi1229032639.ps tmp/84lfi1229032639.png")
> system("convert tmp/942061229032639.ps tmp/942061229032639.png")
> system("convert tmp/10reu01229032639.ps tmp/10reu01229032639.png")
>
>
> proc.time()
user system elapsed
6.021 3.049 6.453