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(114.08,136.49,112.95,142.62,135.31,141.71,134.31,149.51,133.03,147.39,140.11,131.96,124.69,136.38,131.68,127.34,150.95,133.85,137.26,125.14,130.51,141.25,143.15,149.32,118.01,120.92,122.56,134.85,147.97,131.93,135.74,134.22,151.62,143.07,154.82,145.37,145.59,134.32,147.12,126.31,175.86,162.21,140.66,124.09,152.69,153.91,154.38,154.34,132.45,138.70,136.44,150.98,153.24,146.39,154.11,178.30,155.93,168.23,142.53,162.52,148.73,158.86,147.73,152.17,166.79,171.01,144.30,171.49,156.07,189.62,161.70,177.46,152.10,179.98,140.45,156.96,155.56,167.89,174.53,194.78,167.16,192.78,159.48,165.06,173.22,196.60,176.13,151.64,180.31,187.02,185.84,210.99,169.43,219.08,195.25,235.68,174.99,241.44,156.42,187.46,182.08,229.57,182.00,208.44,153.28,215.09,136.72,217.00,130.19,171.08,132.04,178.41,143.89,196.34,133.38,172.11,127.98,154.93,150.45,182.26,133.55,181.74),dim=c(2,61),dimnames=list(c('InvoerEU','InvoerAM'),1:61))
> y <- array(NA,dim=c(2,61),dimnames=list(c('InvoerEU','InvoerAM'),1:61))
> 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
InvoerEU InvoerAM M1 M2 M3 M4 M5 M6 M7 M8 M9 M10 M11
1 114.08 136.49 1 0 0 0 0 0 0 0 0 0 0
2 112.95 142.62 0 1 0 0 0 0 0 0 0 0 0
3 135.31 141.71 0 0 1 0 0 0 0 0 0 0 0
4 134.31 149.51 0 0 0 1 0 0 0 0 0 0 0
5 133.03 147.39 0 0 0 0 1 0 0 0 0 0 0
6 140.11 131.96 0 0 0 0 0 1 0 0 0 0 0
7 124.69 136.38 0 0 0 0 0 0 1 0 0 0 0
8 131.68 127.34 0 0 0 0 0 0 0 1 0 0 0
9 150.95 133.85 0 0 0 0 0 0 0 0 1 0 0
10 137.26 125.14 0 0 0 0 0 0 0 0 0 1 0
11 130.51 141.25 0 0 0 0 0 0 0 0 0 0 1
12 143.15 149.32 0 0 0 0 0 0 0 0 0 0 0
13 118.01 120.92 1 0 0 0 0 0 0 0 0 0 0
14 122.56 134.85 0 1 0 0 0 0 0 0 0 0 0
15 147.97 131.93 0 0 1 0 0 0 0 0 0 0 0
16 135.74 134.22 0 0 0 1 0 0 0 0 0 0 0
17 151.62 143.07 0 0 0 0 1 0 0 0 0 0 0
18 154.82 145.37 0 0 0 0 0 1 0 0 0 0 0
19 145.59 134.32 0 0 0 0 0 0 1 0 0 0 0
20 147.12 126.31 0 0 0 0 0 0 0 1 0 0 0
21 175.86 162.21 0 0 0 0 0 0 0 0 1 0 0
22 140.66 124.09 0 0 0 0 0 0 0 0 0 1 0
23 152.69 153.91 0 0 0 0 0 0 0 0 0 0 1
24 154.38 154.34 0 0 0 0 0 0 0 0 0 0 0
25 132.45 138.70 1 0 0 0 0 0 0 0 0 0 0
26 136.44 150.98 0 1 0 0 0 0 0 0 0 0 0
27 153.24 146.39 0 0 1 0 0 0 0 0 0 0 0
28 154.11 178.30 0 0 0 1 0 0 0 0 0 0 0
29 155.93 168.23 0 0 0 0 1 0 0 0 0 0 0
30 142.53 162.52 0 0 0 0 0 1 0 0 0 0 0
31 148.73 158.86 0 0 0 0 0 0 1 0 0 0 0
32 147.73 152.17 0 0 0 0 0 0 0 1 0 0 0
33 166.79 171.01 0 0 0 0 0 0 0 0 1 0 0
34 144.30 171.49 0 0 0 0 0 0 0 0 0 1 0
35 156.07 189.62 0 0 0 0 0 0 0 0 0 0 1
36 161.70 177.46 0 0 0 0 0 0 0 0 0 0 0
37 152.10 179.98 1 0 0 0 0 0 0 0 0 0 0
38 140.45 156.96 0 1 0 0 0 0 0 0 0 0 0
39 155.56 167.89 0 0 1 0 0 0 0 0 0 0 0
40 174.53 194.78 0 0 0 1 0 0 0 0 0 0 0
41 167.16 192.78 0 0 0 0 1 0 0 0 0 0 0
42 159.48 165.06 0 0 0 0 0 1 0 0 0 0 0
43 173.22 196.60 0 0 0 0 0 0 1 0 0 0 0
44 176.13 151.64 0 0 0 0 0 0 0 1 0 0 0
45 180.31 187.02 0 0 0 0 0 0 0 0 1 0 0
46 185.84 210.99 0 0 0 0 0 0 0 0 0 1 0
47 169.43 219.08 0 0 0 0 0 0 0 0 0 0 1
48 195.25 235.68 0 0 0 0 0 0 0 0 0 0 0
49 174.99 241.44 1 0 0 0 0 0 0 0 0 0 0
50 156.42 187.46 0 1 0 0 0 0 0 0 0 0 0
51 182.08 229.57 0 0 1 0 0 0 0 0 0 0 0
52 182.00 208.44 0 0 0 1 0 0 0 0 0 0 0
53 153.28 215.09 0 0 0 0 1 0 0 0 0 0 0
54 136.72 217.00 0 0 0 0 0 1 0 0 0 0 0
55 130.19 171.08 0 0 0 0 0 0 1 0 0 0 0
56 132.04 178.41 0 0 0 0 0 0 0 1 0 0 0
57 143.89 196.34 0 0 0 0 0 0 0 0 1 0 0
58 133.38 172.11 0 0 0 0 0 0 0 0 0 1 0
59 127.98 154.93 0 0 0 0 0 0 0 0 0 0 1
60 150.45 182.26 0 0 0 0 0 0 0 0 0 0 0
61 133.55 181.74 1 0 0 0 0 0 0 0 0 0 0
> k <- length(x[1,])
> df <- as.data.frame(x)
> (mylm <- lm(df))
Call:
lm(formula = df)
Coefficients:
(Intercept) InvoerAM M1 M2 M3 M4
90.7173 0.3908 -18.2714 -17.3592 0.2213 -2.2055
M5 M6 M7 M8 M9 M10
-6.2419 -8.2241 -8.5440 -1.2914 6.3748 -5.2542
M11
-10.5026
> (mysum <- summary(mylm))
Call:
lm(formula = df)
Residuals:
Min 1Q Median 3Q Max
-30.575 -8.180 2.869 8.333 27.445
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 90.7173 12.2525 7.404 1.77e-09 ***
InvoerAM 0.3908 0.0595 6.568 3.38e-08 ***
M1 -18.2714 8.1239 -2.249 0.0291 *
M2 -17.3592 8.5774 -2.024 0.0486 *
M3 0.2213 8.5006 0.026 0.9793
M4 -2.2055 8.4545 -0.261 0.7953
M5 -6.2419 8.4538 -0.738 0.4639
M6 -8.2241 8.4947 -0.968 0.3378
M7 -8.5440 8.5314 -1.001 0.3216
M8 -1.2914 8.6654 -0.149 0.8822
M9 6.3748 8.4648 0.753 0.4551
M10 -5.2542 8.5207 -0.617 0.5404
M11 -10.5026 8.4585 -1.242 0.2204
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 13.35 on 48 degrees of freedom
Multiple R-squared: 0.5849, Adjusted R-squared: 0.4811
F-statistic: 5.636 on 12 and 48 DF, p-value: 6.456e-06
> 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,] 4.863387e-02 9.726775e-02 0.9513661
[2,] 8.550241e-02 1.710048e-01 0.9144976
[3,] 1.363818e-01 2.727635e-01 0.8636182
[4,] 1.682843e-01 3.365685e-01 0.8317157
[5,] 1.419976e-01 2.839951e-01 0.8580024
[6,] 1.940634e-01 3.881268e-01 0.8059366
[7,] 1.256861e-01 2.513721e-01 0.8743139
[8,] 1.215595e-01 2.431190e-01 0.8784405
[9,] 8.183330e-02 1.636666e-01 0.9181667
[10,] 6.108483e-02 1.221697e-01 0.9389152
[11,] 4.505132e-02 9.010263e-02 0.9549487
[12,] 2.688590e-02 5.377180e-02 0.9731141
[13,] 1.646082e-02 3.292165e-02 0.9835392
[14,] 8.891208e-03 1.778242e-02 0.9911088
[15,] 9.224899e-03 1.844980e-02 0.9907751
[16,] 4.763358e-03 9.526716e-03 0.9952366
[17,] 2.321914e-03 4.643828e-03 0.9976781
[18,] 1.363068e-03 2.726135e-03 0.9986369
[19,] 9.578224e-04 1.915645e-03 0.9990422
[20,] 4.151666e-04 8.303333e-04 0.9995848
[21,] 1.714720e-04 3.429439e-04 0.9998285
[22,] 1.145675e-04 2.291350e-04 0.9998854
[23,] 5.849249e-05 1.169850e-04 0.9999415
[24,] 2.022154e-05 4.044307e-05 0.9999798
[25,] 1.306287e-05 2.612575e-05 0.9999869
[26,] 6.765693e-06 1.353139e-05 0.9999932
[27,] 3.980261e-05 7.960521e-05 0.9999602
[28,] 3.307203e-05 6.614406e-05 0.9999669
[29,] 1.591329e-02 3.182657e-02 0.9840867
[30,] 3.316432e-01 6.632864e-01 0.6683568
> postscript(file="/var/www/html/rcomp/tmp/1w54d1259086654.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/2ytj21259086654.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/39at61259086654.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/4849v1259086654.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/51ohs1259086654.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 = 61
Frequency = 1
1 2 3 4 5 6
-11.7048072 -16.1424966 -11.0074679 -12.6288039 -9.0439422 6.0481927
7 8 9 10 11 12
-10.7792568 -7.5090716 1.5506657 2.8935022 -4.9037793 -5.9200320
13 14 15 16 17 18
-1.6902072 -3.4960584 5.4744581 -5.2236251 11.2342706 15.5176991
19 20 21 22 23 24
10.9257706 8.3334421 15.3778619 6.7038317 12.3288195 3.3482024
25 26 27 28 29 30
5.8015470 4.0804992 5.0936350 -4.0796474 5.7119947 -3.4743491
31 32 33 34 35 36
4.4757845 -1.1623867 2.8689101 -8.1796133 1.7537094 1.6331380
37 38 39 40 41 42
9.3197366 5.7535751 -0.9883496 9.9001338 7.3481007 12.4830444
43 44 45 46 47 48
14.2173707 27.4447320 10.1323625 17.9241825 3.6010366 12.4313454
49 50 51 52 53 54
8.1917845 9.8044807 1.4277245 12.0319426 -15.2504238 -30.5745871
55 56 57 58 59 60
-18.8396690 -27.1067158 -29.9298001 -19.3419031 -12.7797862 -11.4926539
61
-9.9180537
> postscript(file="/var/www/html/rcomp/tmp/6xzmo1259086654.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 = 61
Frequency = 1
lag(myerror, k = 1) myerror
0 -11.7048072 NA
1 -16.1424966 -11.7048072
2 -11.0074679 -16.1424966
3 -12.6288039 -11.0074679
4 -9.0439422 -12.6288039
5 6.0481927 -9.0439422
6 -10.7792568 6.0481927
7 -7.5090716 -10.7792568
8 1.5506657 -7.5090716
9 2.8935022 1.5506657
10 -4.9037793 2.8935022
11 -5.9200320 -4.9037793
12 -1.6902072 -5.9200320
13 -3.4960584 -1.6902072
14 5.4744581 -3.4960584
15 -5.2236251 5.4744581
16 11.2342706 -5.2236251
17 15.5176991 11.2342706
18 10.9257706 15.5176991
19 8.3334421 10.9257706
20 15.3778619 8.3334421
21 6.7038317 15.3778619
22 12.3288195 6.7038317
23 3.3482024 12.3288195
24 5.8015470 3.3482024
25 4.0804992 5.8015470
26 5.0936350 4.0804992
27 -4.0796474 5.0936350
28 5.7119947 -4.0796474
29 -3.4743491 5.7119947
30 4.4757845 -3.4743491
31 -1.1623867 4.4757845
32 2.8689101 -1.1623867
33 -8.1796133 2.8689101
34 1.7537094 -8.1796133
35 1.6331380 1.7537094
36 9.3197366 1.6331380
37 5.7535751 9.3197366
38 -0.9883496 5.7535751
39 9.9001338 -0.9883496
40 7.3481007 9.9001338
41 12.4830444 7.3481007
42 14.2173707 12.4830444
43 27.4447320 14.2173707
44 10.1323625 27.4447320
45 17.9241825 10.1323625
46 3.6010366 17.9241825
47 12.4313454 3.6010366
48 8.1917845 12.4313454
49 9.8044807 8.1917845
50 1.4277245 9.8044807
51 12.0319426 1.4277245
52 -15.2504238 12.0319426
53 -30.5745871 -15.2504238
54 -18.8396690 -30.5745871
55 -27.1067158 -18.8396690
56 -29.9298001 -27.1067158
57 -19.3419031 -29.9298001
58 -12.7797862 -19.3419031
59 -11.4926539 -12.7797862
60 -9.9180537 -11.4926539
61 NA -9.9180537
> dum1 <- dum[2:length(myerror),]
> dum1
lag(myerror, k = 1) myerror
[1,] -16.1424966 -11.7048072
[2,] -11.0074679 -16.1424966
[3,] -12.6288039 -11.0074679
[4,] -9.0439422 -12.6288039
[5,] 6.0481927 -9.0439422
[6,] -10.7792568 6.0481927
[7,] -7.5090716 -10.7792568
[8,] 1.5506657 -7.5090716
[9,] 2.8935022 1.5506657
[10,] -4.9037793 2.8935022
[11,] -5.9200320 -4.9037793
[12,] -1.6902072 -5.9200320
[13,] -3.4960584 -1.6902072
[14,] 5.4744581 -3.4960584
[15,] -5.2236251 5.4744581
[16,] 11.2342706 -5.2236251
[17,] 15.5176991 11.2342706
[18,] 10.9257706 15.5176991
[19,] 8.3334421 10.9257706
[20,] 15.3778619 8.3334421
[21,] 6.7038317 15.3778619
[22,] 12.3288195 6.7038317
[23,] 3.3482024 12.3288195
[24,] 5.8015470 3.3482024
[25,] 4.0804992 5.8015470
[26,] 5.0936350 4.0804992
[27,] -4.0796474 5.0936350
[28,] 5.7119947 -4.0796474
[29,] -3.4743491 5.7119947
[30,] 4.4757845 -3.4743491
[31,] -1.1623867 4.4757845
[32,] 2.8689101 -1.1623867
[33,] -8.1796133 2.8689101
[34,] 1.7537094 -8.1796133
[35,] 1.6331380 1.7537094
[36,] 9.3197366 1.6331380
[37,] 5.7535751 9.3197366
[38,] -0.9883496 5.7535751
[39,] 9.9001338 -0.9883496
[40,] 7.3481007 9.9001338
[41,] 12.4830444 7.3481007
[42,] 14.2173707 12.4830444
[43,] 27.4447320 14.2173707
[44,] 10.1323625 27.4447320
[45,] 17.9241825 10.1323625
[46,] 3.6010366 17.9241825
[47,] 12.4313454 3.6010366
[48,] 8.1917845 12.4313454
[49,] 9.8044807 8.1917845
[50,] 1.4277245 9.8044807
[51,] 12.0319426 1.4277245
[52,] -15.2504238 12.0319426
[53,] -30.5745871 -15.2504238
[54,] -18.8396690 -30.5745871
[55,] -27.1067158 -18.8396690
[56,] -29.9298001 -27.1067158
[57,] -19.3419031 -29.9298001
[58,] -12.7797862 -19.3419031
[59,] -11.4926539 -12.7797862
[60,] -9.9180537 -11.4926539
> z <- as.data.frame(dum1)
> z
lag(myerror, k = 1) myerror
1 -16.1424966 -11.7048072
2 -11.0074679 -16.1424966
3 -12.6288039 -11.0074679
4 -9.0439422 -12.6288039
5 6.0481927 -9.0439422
6 -10.7792568 6.0481927
7 -7.5090716 -10.7792568
8 1.5506657 -7.5090716
9 2.8935022 1.5506657
10 -4.9037793 2.8935022
11 -5.9200320 -4.9037793
12 -1.6902072 -5.9200320
13 -3.4960584 -1.6902072
14 5.4744581 -3.4960584
15 -5.2236251 5.4744581
16 11.2342706 -5.2236251
17 15.5176991 11.2342706
18 10.9257706 15.5176991
19 8.3334421 10.9257706
20 15.3778619 8.3334421
21 6.7038317 15.3778619
22 12.3288195 6.7038317
23 3.3482024 12.3288195
24 5.8015470 3.3482024
25 4.0804992 5.8015470
26 5.0936350 4.0804992
27 -4.0796474 5.0936350
28 5.7119947 -4.0796474
29 -3.4743491 5.7119947
30 4.4757845 -3.4743491
31 -1.1623867 4.4757845
32 2.8689101 -1.1623867
33 -8.1796133 2.8689101
34 1.7537094 -8.1796133
35 1.6331380 1.7537094
36 9.3197366 1.6331380
37 5.7535751 9.3197366
38 -0.9883496 5.7535751
39 9.9001338 -0.9883496
40 7.3481007 9.9001338
41 12.4830444 7.3481007
42 14.2173707 12.4830444
43 27.4447320 14.2173707
44 10.1323625 27.4447320
45 17.9241825 10.1323625
46 3.6010366 17.9241825
47 12.4313454 3.6010366
48 8.1917845 12.4313454
49 9.8044807 8.1917845
50 1.4277245 9.8044807
51 12.0319426 1.4277245
52 -15.2504238 12.0319426
53 -30.5745871 -15.2504238
54 -18.8396690 -30.5745871
55 -27.1067158 -18.8396690
56 -29.9298001 -27.1067158
57 -19.3419031 -29.9298001
58 -12.7797862 -19.3419031
59 -11.4926539 -12.7797862
60 -9.9180537 -11.4926539
> 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/7s6q91259086654.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/800pl1259086654.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/9i0i01259086654.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/10knxx1259086654.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/11n0we1259086654.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/12o4gz1259086655.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/13f8zb1259086655.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/14htet1259086655.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/15gshg1259086655.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/16nv081259086655.tab")
+ }
>
> system("convert tmp/1w54d1259086654.ps tmp/1w54d1259086654.png")
> system("convert tmp/2ytj21259086654.ps tmp/2ytj21259086654.png")
> system("convert tmp/39at61259086654.ps tmp/39at61259086654.png")
> system("convert tmp/4849v1259086654.ps tmp/4849v1259086654.png")
> system("convert tmp/51ohs1259086654.ps tmp/51ohs1259086654.png")
> system("convert tmp/6xzmo1259086654.ps tmp/6xzmo1259086654.png")
> system("convert tmp/7s6q91259086654.ps tmp/7s6q91259086654.png")
> system("convert tmp/800pl1259086654.ps tmp/800pl1259086654.png")
> system("convert tmp/9i0i01259086654.ps tmp/9i0i01259086654.png")
> system("convert tmp/10knxx1259086654.ps tmp/10knxx1259086654.png")
>
>
> proc.time()
user system elapsed
2.461 1.593 3.519