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(10.9,96.8,10,114.1,9.2,110.3,9.2,103.9,9.5,101.6,9.6,94.6,9.5,95.9,9.1,104.7,8.9,102.8,9,98.1,10.1,113.9,10.3,80.9,10.2,95.7,9.6,113.2,9.2,105.9,9.3,108.8,9.4,102.3,9.4,99,9.2,100.7,9,115.5,9,100.7,9,109.9,9.8,114.6,10,85.4,9.8,100.5,9.3,114.8,9,116.5,9,112.9,9.1,102,9.1,106,9.1,105.3,9.2,118.8,8.8,106.1,8.3,109.3,8.4,117.2,8.1,92.5,7.7,104.2,7.9,112.5,7.9,122.4,8,113.3,7.9,100,7.6,110.7,7.1,112.8,6.8,109.8,6.5,117.3,6.9,109.1,8.2,115.9,8.7,96,8.3,99.8,7.9,116.8,7.5,115.7,7.8,99.4,8.3,94.3,8.4,91,8.2,93.2,7.7,103.1,7.2,94.1,7.3,91.8,8.1,102.7,8.5,82.6),dim=c(2,60),dimnames=list(c('Y','X'),1:60))
> y <- array(NA,dim=c(2,60),dimnames=list(c('Y','X'),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 = '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
Y X M1 M2 M3 M4 M5 M6 M7 M8 M9 M10 M11
1 10.9 96.8 1 0 0 0 0 0 0 0 0 0 0
2 10.0 114.1 0 1 0 0 0 0 0 0 0 0 0
3 9.2 110.3 0 0 1 0 0 0 0 0 0 0 0
4 9.2 103.9 0 0 0 1 0 0 0 0 0 0 0
5 9.5 101.6 0 0 0 0 1 0 0 0 0 0 0
6 9.6 94.6 0 0 0 0 0 1 0 0 0 0 0
7 9.5 95.9 0 0 0 0 0 0 1 0 0 0 0
8 9.1 104.7 0 0 0 0 0 0 0 1 0 0 0
9 8.9 102.8 0 0 0 0 0 0 0 0 1 0 0
10 9.0 98.1 0 0 0 0 0 0 0 0 0 1 0
11 10.1 113.9 0 0 0 0 0 0 0 0 0 0 1
12 10.3 80.9 0 0 0 0 0 0 0 0 0 0 0
13 10.2 95.7 1 0 0 0 0 0 0 0 0 0 0
14 9.6 113.2 0 1 0 0 0 0 0 0 0 0 0
15 9.2 105.9 0 0 1 0 0 0 0 0 0 0 0
16 9.3 108.8 0 0 0 1 0 0 0 0 0 0 0
17 9.4 102.3 0 0 0 0 1 0 0 0 0 0 0
18 9.4 99.0 0 0 0 0 0 1 0 0 0 0 0
19 9.2 100.7 0 0 0 0 0 0 1 0 0 0 0
20 9.0 115.5 0 0 0 0 0 0 0 1 0 0 0
21 9.0 100.7 0 0 0 0 0 0 0 0 1 0 0
22 9.0 109.9 0 0 0 0 0 0 0 0 0 1 0
23 9.8 114.6 0 0 0 0 0 0 0 0 0 0 1
24 10.0 85.4 0 0 0 0 0 0 0 0 0 0 0
25 9.8 100.5 1 0 0 0 0 0 0 0 0 0 0
26 9.3 114.8 0 1 0 0 0 0 0 0 0 0 0
27 9.0 116.5 0 0 1 0 0 0 0 0 0 0 0
28 9.0 112.9 0 0 0 1 0 0 0 0 0 0 0
29 9.1 102.0 0 0 0 0 1 0 0 0 0 0 0
30 9.1 106.0 0 0 0 0 0 1 0 0 0 0 0
31 9.1 105.3 0 0 0 0 0 0 1 0 0 0 0
32 9.2 118.8 0 0 0 0 0 0 0 1 0 0 0
33 8.8 106.1 0 0 0 0 0 0 0 0 1 0 0
34 8.3 109.3 0 0 0 0 0 0 0 0 0 1 0
35 8.4 117.2 0 0 0 0 0 0 0 0 0 0 1
36 8.1 92.5 0 0 0 0 0 0 0 0 0 0 0
37 7.7 104.2 1 0 0 0 0 0 0 0 0 0 0
38 7.9 112.5 0 1 0 0 0 0 0 0 0 0 0
39 7.9 122.4 0 0 1 0 0 0 0 0 0 0 0
40 8.0 113.3 0 0 0 1 0 0 0 0 0 0 0
41 7.9 100.0 0 0 0 0 1 0 0 0 0 0 0
42 7.6 110.7 0 0 0 0 0 1 0 0 0 0 0
43 7.1 112.8 0 0 0 0 0 0 1 0 0 0 0
44 6.8 109.8 0 0 0 0 0 0 0 1 0 0 0
45 6.5 117.3 0 0 0 0 0 0 0 0 1 0 0
46 6.9 109.1 0 0 0 0 0 0 0 0 0 1 0
47 8.2 115.9 0 0 0 0 0 0 0 0 0 0 1
48 8.7 96.0 0 0 0 0 0 0 0 0 0 0 0
49 8.3 99.8 1 0 0 0 0 0 0 0 0 0 0
50 7.9 116.8 0 1 0 0 0 0 0 0 0 0 0
51 7.5 115.7 0 0 1 0 0 0 0 0 0 0 0
52 7.8 99.4 0 0 0 1 0 0 0 0 0 0 0
53 8.3 94.3 0 0 0 0 1 0 0 0 0 0 0
54 8.4 91.0 0 0 0 0 0 1 0 0 0 0 0
55 8.2 93.2 0 0 0 0 0 0 1 0 0 0 0
56 7.7 103.1 0 0 0 0 0 0 0 1 0 0 0
57 7.2 94.1 0 0 0 0 0 0 0 0 1 0 0
58 7.3 91.8 0 0 0 0 0 0 0 0 0 1 0
59 8.1 102.7 0 0 0 0 0 0 0 0 0 0 1
60 8.5 82.6 0 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) X M1 M2 M3 M4
11.66897 -0.02914 0.60732 0.60089 0.21739 0.12800
M5 M6 M7 M8 M9 M10
0.08597 0.07238 -0.08916 -0.09275 -0.55282 -0.54913
M11
0.53951
> (mysum <- summary(mylm))
Call:
lm(formula = df)
Residuals:
Min 1Q Median 3Q Max
-1.5769 -0.8842 0.4149 0.7071 1.4442
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 11.66897 1.94794 5.990 2.78e-07 ***
X -0.02914 0.02172 -1.342 0.186
M1 0.60732 0.66051 0.919 0.363
M2 0.60089 0.84145 0.714 0.479
M3 0.21739 0.83965 0.259 0.797
M4 0.12800 0.74923 0.171 0.865
M5 0.08597 0.66608 0.129 0.898
M6 0.07238 0.66805 0.108 0.914
M7 -0.08916 0.68046 -0.131 0.896
M8 -0.09275 0.78525 -0.118 0.906
M9 -0.55282 0.70790 -0.781 0.439
M10 -0.54913 0.70173 -0.783 0.438
M11 0.53951 0.82042 0.658 0.514
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 0.9608 on 47 degrees of freedom
Multiple R-squared: 0.1872, Adjusted R-squared: -0.02031
F-statistic: 0.9022 on 12 and 47 DF, p-value: 0.5515
> 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.983025e-02 9.966050e-02 0.95016975
[2,] 1.464181e-02 2.928361e-02 0.98535819
[3,] 5.223029e-03 1.044606e-02 0.99477697
[4,] 1.939670e-03 3.879341e-03 0.99806033
[5,] 5.191759e-04 1.038352e-03 0.99948082
[6,] 1.587658e-04 3.175316e-04 0.99984123
[7,] 4.981656e-05 9.963313e-05 0.99995018
[8,] 2.657940e-05 5.315880e-05 0.99997342
[9,] 1.452437e-05 2.904874e-05 0.99998548
[10,] 1.058118e-04 2.116237e-04 0.99989419
[11,] 1.380512e-04 2.761024e-04 0.99986195
[12,] 7.369070e-05 1.473814e-04 0.99992631
[13,] 3.599865e-05 7.199729e-05 0.99996400
[14,] 2.814324e-05 5.628648e-05 0.99997186
[15,] 1.711292e-05 3.422583e-05 0.99998289
[16,] 1.684719e-05 3.369437e-05 0.99998315
[17,] 2.318622e-04 4.637245e-04 0.99976814
[18,] 5.987369e-03 1.197474e-02 0.99401263
[19,] 1.253684e-01 2.507369e-01 0.87463157
[20,] 5.834632e-01 8.330736e-01 0.41653682
[21,] 8.437609e-01 3.124783e-01 0.15623914
[22,] 9.662026e-01 6.759481e-02 0.03379741
[23,] 9.731224e-01 5.375511e-02 0.02687756
[24,] 9.675237e-01 6.495259e-02 0.03247630
[25,] 9.672732e-01 6.545360e-02 0.03272680
[26,] 9.545405e-01 9.091890e-02 0.04545945
[27,] 9.142840e-01 1.714319e-01 0.08571596
[28,] 9.022732e-01 1.954536e-01 0.09772679
[29,] 9.574172e-01 8.516558e-02 0.04258279
> postscript(file="/var/www/html/rcomp/tmp/1mf1b1258710219.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/2x97t1258710219.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/3c0391258710219.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/46xpz1258710219.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/5cach1258710219.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 7
1.4442420 1.0547552 0.5275285 0.4304423 0.7054548 0.6150807 0.7144979
8 9 10 11 12 13 14
0.5744979 0.7792072 0.7385772 1.2103032 0.9882740 0.7121905 0.6285313
15 16 17 18 19 20 21
0.3993227 0.6732170 0.6258512 0.5432865 0.5543588 0.7891850 0.8180181
22 23 24 25 26 27 28
1.0824020 0.9306996 0.8193936 0.4520515 0.3751516 0.5081822 0.4926815
29 30 31 32 33 34 35
0.3171099 0.4472504 0.5883922 1.0853394 0.7753616 0.3649193 -0.3935424
36 37 38 39 40 41 42
-0.8737288 -1.5401391 -1.0918651 -0.4199054 -0.4956634 -0.9411655 -0.9158025
43 44 45 46 47 48 49
-1.1930751 -1.5768999 -1.1982962 -1.0409082 -0.6314214 -0.1717469 -1.0683449
50 51 52 53 54 55 56
-0.9665730 -1.0151280 -1.1006773 -0.7072504 -0.6898150 -0.6641739 -0.8721224
57 58 59 60
-1.1742907 -1.1449903 -1.1160390 -0.7621919
> postscript(file="/var/www/html/rcomp/tmp/6ad031258710220.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 1.4442420 NA
1 1.0547552 1.4442420
2 0.5275285 1.0547552
3 0.4304423 0.5275285
4 0.7054548 0.4304423
5 0.6150807 0.7054548
6 0.7144979 0.6150807
7 0.5744979 0.7144979
8 0.7792072 0.5744979
9 0.7385772 0.7792072
10 1.2103032 0.7385772
11 0.9882740 1.2103032
12 0.7121905 0.9882740
13 0.6285313 0.7121905
14 0.3993227 0.6285313
15 0.6732170 0.3993227
16 0.6258512 0.6732170
17 0.5432865 0.6258512
18 0.5543588 0.5432865
19 0.7891850 0.5543588
20 0.8180181 0.7891850
21 1.0824020 0.8180181
22 0.9306996 1.0824020
23 0.8193936 0.9306996
24 0.4520515 0.8193936
25 0.3751516 0.4520515
26 0.5081822 0.3751516
27 0.4926815 0.5081822
28 0.3171099 0.4926815
29 0.4472504 0.3171099
30 0.5883922 0.4472504
31 1.0853394 0.5883922
32 0.7753616 1.0853394
33 0.3649193 0.7753616
34 -0.3935424 0.3649193
35 -0.8737288 -0.3935424
36 -1.5401391 -0.8737288
37 -1.0918651 -1.5401391
38 -0.4199054 -1.0918651
39 -0.4956634 -0.4199054
40 -0.9411655 -0.4956634
41 -0.9158025 -0.9411655
42 -1.1930751 -0.9158025
43 -1.5768999 -1.1930751
44 -1.1982962 -1.5768999
45 -1.0409082 -1.1982962
46 -0.6314214 -1.0409082
47 -0.1717469 -0.6314214
48 -1.0683449 -0.1717469
49 -0.9665730 -1.0683449
50 -1.0151280 -0.9665730
51 -1.1006773 -1.0151280
52 -0.7072504 -1.1006773
53 -0.6898150 -0.7072504
54 -0.6641739 -0.6898150
55 -0.8721224 -0.6641739
56 -1.1742907 -0.8721224
57 -1.1449903 -1.1742907
58 -1.1160390 -1.1449903
59 -0.7621919 -1.1160390
60 NA -0.7621919
> dum1 <- dum[2:length(myerror),]
> dum1
lag(myerror, k = 1) myerror
[1,] 1.0547552 1.4442420
[2,] 0.5275285 1.0547552
[3,] 0.4304423 0.5275285
[4,] 0.7054548 0.4304423
[5,] 0.6150807 0.7054548
[6,] 0.7144979 0.6150807
[7,] 0.5744979 0.7144979
[8,] 0.7792072 0.5744979
[9,] 0.7385772 0.7792072
[10,] 1.2103032 0.7385772
[11,] 0.9882740 1.2103032
[12,] 0.7121905 0.9882740
[13,] 0.6285313 0.7121905
[14,] 0.3993227 0.6285313
[15,] 0.6732170 0.3993227
[16,] 0.6258512 0.6732170
[17,] 0.5432865 0.6258512
[18,] 0.5543588 0.5432865
[19,] 0.7891850 0.5543588
[20,] 0.8180181 0.7891850
[21,] 1.0824020 0.8180181
[22,] 0.9306996 1.0824020
[23,] 0.8193936 0.9306996
[24,] 0.4520515 0.8193936
[25,] 0.3751516 0.4520515
[26,] 0.5081822 0.3751516
[27,] 0.4926815 0.5081822
[28,] 0.3171099 0.4926815
[29,] 0.4472504 0.3171099
[30,] 0.5883922 0.4472504
[31,] 1.0853394 0.5883922
[32,] 0.7753616 1.0853394
[33,] 0.3649193 0.7753616
[34,] -0.3935424 0.3649193
[35,] -0.8737288 -0.3935424
[36,] -1.5401391 -0.8737288
[37,] -1.0918651 -1.5401391
[38,] -0.4199054 -1.0918651
[39,] -0.4956634 -0.4199054
[40,] -0.9411655 -0.4956634
[41,] -0.9158025 -0.9411655
[42,] -1.1930751 -0.9158025
[43,] -1.5768999 -1.1930751
[44,] -1.1982962 -1.5768999
[45,] -1.0409082 -1.1982962
[46,] -0.6314214 -1.0409082
[47,] -0.1717469 -0.6314214
[48,] -1.0683449 -0.1717469
[49,] -0.9665730 -1.0683449
[50,] -1.0151280 -0.9665730
[51,] -1.1006773 -1.0151280
[52,] -0.7072504 -1.1006773
[53,] -0.6898150 -0.7072504
[54,] -0.6641739 -0.6898150
[55,] -0.8721224 -0.6641739
[56,] -1.1742907 -0.8721224
[57,] -1.1449903 -1.1742907
[58,] -1.1160390 -1.1449903
[59,] -0.7621919 -1.1160390
> z <- as.data.frame(dum1)
> z
lag(myerror, k = 1) myerror
1 1.0547552 1.4442420
2 0.5275285 1.0547552
3 0.4304423 0.5275285
4 0.7054548 0.4304423
5 0.6150807 0.7054548
6 0.7144979 0.6150807
7 0.5744979 0.7144979
8 0.7792072 0.5744979
9 0.7385772 0.7792072
10 1.2103032 0.7385772
11 0.9882740 1.2103032
12 0.7121905 0.9882740
13 0.6285313 0.7121905
14 0.3993227 0.6285313
15 0.6732170 0.3993227
16 0.6258512 0.6732170
17 0.5432865 0.6258512
18 0.5543588 0.5432865
19 0.7891850 0.5543588
20 0.8180181 0.7891850
21 1.0824020 0.8180181
22 0.9306996 1.0824020
23 0.8193936 0.9306996
24 0.4520515 0.8193936
25 0.3751516 0.4520515
26 0.5081822 0.3751516
27 0.4926815 0.5081822
28 0.3171099 0.4926815
29 0.4472504 0.3171099
30 0.5883922 0.4472504
31 1.0853394 0.5883922
32 0.7753616 1.0853394
33 0.3649193 0.7753616
34 -0.3935424 0.3649193
35 -0.8737288 -0.3935424
36 -1.5401391 -0.8737288
37 -1.0918651 -1.5401391
38 -0.4199054 -1.0918651
39 -0.4956634 -0.4199054
40 -0.9411655 -0.4956634
41 -0.9158025 -0.9411655
42 -1.1930751 -0.9158025
43 -1.5768999 -1.1930751
44 -1.1982962 -1.5768999
45 -1.0409082 -1.1982962
46 -0.6314214 -1.0409082
47 -0.1717469 -0.6314214
48 -1.0683449 -0.1717469
49 -0.9665730 -1.0683449
50 -1.0151280 -0.9665730
51 -1.1006773 -1.0151280
52 -0.7072504 -1.1006773
53 -0.6898150 -0.7072504
54 -0.6641739 -0.6898150
55 -0.8721224 -0.6641739
56 -1.1742907 -0.8721224
57 -1.1449903 -1.1742907
58 -1.1160390 -1.1449903
59 -0.7621919 -1.1160390
> 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/7jpzm1258710220.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/8m7fh1258710220.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/9d4th1258710220.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/10x22r1258710220.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/11mzhk1258710220.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/12aas51258710220.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/13pgq51258710220.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/14vw9p1258710220.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/15u9rp1258710220.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/16naof1258710220.tab")
+ }
>
> system("convert tmp/1mf1b1258710219.ps tmp/1mf1b1258710219.png")
> system("convert tmp/2x97t1258710219.ps tmp/2x97t1258710219.png")
> system("convert tmp/3c0391258710219.ps tmp/3c0391258710219.png")
> system("convert tmp/46xpz1258710219.ps tmp/46xpz1258710219.png")
> system("convert tmp/5cach1258710219.ps tmp/5cach1258710219.png")
> system("convert tmp/6ad031258710220.ps tmp/6ad031258710220.png")
> system("convert tmp/7jpzm1258710220.ps tmp/7jpzm1258710220.png")
> system("convert tmp/8m7fh1258710220.ps tmp/8m7fh1258710220.png")
> system("convert tmp/9d4th1258710220.ps tmp/9d4th1258710220.png")
> system("convert tmp/10x22r1258710220.ps tmp/10x22r1258710220.png")
>
>
> proc.time()
user system elapsed
2.380 1.533 3.212