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(1469798.00
+ ,10467.48
+ ,1368839.00
+ ,1207763.00
+ ,1008380.00
+ ,989236.00
+ ,1498721.00
+ ,10274.97
+ ,1469798.00
+ ,1368839.00
+ ,1207763.00
+ ,1008380.00
+ ,1761769.00
+ ,10640.91
+ ,1498721.00
+ ,1469798.00
+ ,1368839.00
+ ,1207763.00
+ ,1653214.00
+ ,10481.60
+ ,1761769.00
+ ,1498721.00
+ ,1469798.00
+ ,1368839.00
+ ,1599104.00
+ ,10568.70
+ ,1653214.00
+ ,1761769.00
+ ,1498721.00
+ ,1469798.00
+ ,1421179.00
+ ,10440.07
+ ,1599104.00
+ ,1653214.00
+ ,1761769.00
+ ,1498721.00
+ ,1163995.00
+ ,10805.87
+ ,1421179.00
+ ,1599104.00
+ ,1653214.00
+ ,1761769.00
+ ,1037735.00
+ ,10717.50
+ ,1163995.00
+ ,1421179.00
+ ,1599104.00
+ ,1653214.00
+ ,1015407.00
+ ,10864.86
+ ,1037735.00
+ ,1163995.00
+ ,1421179.00
+ ,1599104.00
+ ,1039210.00
+ ,10993.41
+ ,1015407.00
+ ,1037735.00
+ ,1163995.00
+ ,1421179.00
+ ,1258049.00
+ ,11109.32
+ ,1039210.00
+ ,1015407.00
+ ,1037735.00
+ ,1163995.00
+ ,1469445.00
+ ,11367.14
+ ,1258049.00
+ ,1039210.00
+ ,1015407.00
+ ,1037735.00
+ ,1552346.00
+ ,11168.31
+ ,1469445.00
+ ,1258049.00
+ ,1039210.00
+ ,1015407.00
+ ,1549144.00
+ ,11150.22
+ ,1552346.00
+ ,1469445.00
+ ,1258049.00
+ ,1039210.00
+ ,1785895.00
+ ,11185.68
+ ,1549144.00
+ ,1552346.00
+ ,1469445.00
+ ,1258049.00
+ ,1662335.00
+ ,11381.15
+ ,1785895.00
+ ,1549144.00
+ ,1552346.00
+ ,1469445.00
+ ,1629440.00
+ ,11679.07
+ ,1662335.00
+ ,1785895.00
+ ,1549144.00
+ ,1552346.00
+ ,1467430.00
+ ,12080.73
+ ,1629440.00
+ ,1662335.00
+ ,1785895.00
+ ,1549144.00
+ ,1202209.00
+ ,12221.93
+ ,1467430.00
+ ,1629440.00
+ ,1662335.00
+ ,1785895.00
+ ,1076982.00
+ ,12463.15
+ ,1202209.00
+ ,1467430.00
+ ,1629440.00
+ ,1662335.00
+ ,1039367.00
+ ,12621.69
+ ,1076982.00
+ ,1202209.00
+ ,1467430.00
+ ,1629440.00
+ ,1063449.00
+ ,12268.63
+ ,1039367.00
+ ,1076982.00
+ ,1202209.00
+ ,1467430.00
+ ,1335135.00
+ ,12354.35
+ ,1063449.00
+ ,1039367.00
+ ,1076982.00
+ ,1202209.00
+ ,1491602.00
+ ,13062.91
+ ,1335135.00
+ ,1063449.00
+ ,1039367.00
+ ,1076982.00
+ ,1591972.00
+ ,13627.64
+ ,1491602.00
+ ,1335135.00
+ ,1063449.00
+ ,1039367.00
+ ,1641248.00
+ ,13408.62
+ ,1591972.00
+ ,1491602.00
+ ,1335135.00
+ ,1063449.00
+ ,1898849.00
+ ,13211.99
+ ,1641248.00
+ ,1591972.00
+ ,1491602.00
+ ,1335135.00
+ ,1798580.00
+ ,13357.74
+ ,1898849.00
+ ,1641248.00
+ ,1591972.00
+ ,1491602.00
+ ,1762444.00
+ ,13895.63
+ ,1798580.00
+ ,1898849.00
+ ,1641248.00
+ ,1591972.00
+ ,1622044.00
+ ,13930.01
+ ,1762444.00
+ ,1798580.00
+ ,1898849.00
+ ,1641248.00
+ ,1368955.00
+ ,13371.72
+ ,1622044.00
+ ,1762444.00
+ ,1798580.00
+ ,1898849.00
+ ,1262973.00
+ ,13264.82
+ ,1368955.00
+ ,1622044.00
+ ,1762444.00
+ ,1798580.00
+ ,1195650.00
+ ,12650.36
+ ,1262973.00
+ ,1368955.00
+ ,1622044.00
+ ,1762444.00
+ ,1269530.00
+ ,12266.39
+ ,1195650.00
+ ,1262973.00
+ ,1368955.00
+ ,1622044.00
+ ,1479279.00
+ ,12262.89
+ ,1269530.00
+ ,1195650.00
+ ,1262973.00
+ ,1368955.00
+ ,1607819.00
+ ,12820.13
+ ,1479279.00
+ ,1269530.00
+ ,1195650.00
+ ,1262973.00
+ ,1712466.00
+ ,12638.32
+ ,1607819.00
+ ,1479279.00
+ ,1269530.00
+ ,1195650.00
+ ,1721766.00
+ ,11350.01
+ ,1712466.00
+ ,1607819.00
+ ,1479279.00
+ ,1269530.00
+ ,1949843.00
+ ,11378.02
+ ,1721766.00
+ ,1712466.00
+ ,1607819.00
+ ,1479279.00
+ ,1821326.00
+ ,11543.55
+ ,1949843.00
+ ,1721766.00
+ ,1712466.00
+ ,1607819.00
+ ,1757802.00
+ ,10850.66
+ ,1821326.00
+ ,1949843.00
+ ,1721766.00
+ ,1712466.00
+ ,1590367.00
+ ,9325.01
+ ,1757802.00
+ ,1821326.00
+ ,1949843.00
+ ,1721766.00
+ ,1260647.00
+ ,8829.04
+ ,1590367.00
+ ,1757802.00
+ ,1821326.00
+ ,1949843.00
+ ,1149235.00
+ ,8776.39
+ ,1260647.00
+ ,1590367.00
+ ,1757802.00
+ ,1821326.00
+ ,1016367.00
+ ,8000.86
+ ,1149235.00
+ ,1260647.00
+ ,1590367.00
+ ,1757802.00
+ ,1027885.00
+ ,7062.93
+ ,1016367.00
+ ,1149235.00
+ ,1260647.00
+ ,1590367.00
+ ,1262159.00
+ ,7608.92
+ ,1027885.00
+ ,1016367.00
+ ,1149235.00
+ ,1260647.00
+ ,1520854.00
+ ,8168.12
+ ,1262159.00
+ ,1027885.00
+ ,1016367.00
+ ,1149235.00
+ ,1544144.00
+ ,8500.33
+ ,1520854.00
+ ,1262159.00
+ ,1027885.00
+ ,1016367.00
+ ,1564709.00
+ ,8447.00
+ ,1544144.00
+ ,1520854.00
+ ,1262159.00
+ ,1027885.00
+ ,1821776.00
+ ,9171.61
+ ,1564709.00
+ ,1544144.00
+ ,1520854.00
+ ,1262159.00
+ ,1741365.00
+ ,9496.28
+ ,1821776.00
+ ,1564709.00
+ ,1544144.00
+ ,1520854.00
+ ,1623386.00
+ ,9712.28
+ ,1741365.00
+ ,1821776.00
+ ,1564709.00
+ ,1544144.00
+ ,1498658.00
+ ,9712.73
+ ,1623386.00
+ ,1741365.00
+ ,1821776.00
+ ,1564709.00
+ ,1241822.00
+ ,10344.84
+ ,1498658.00
+ ,1623386.00
+ ,1741365.00
+ ,1821776.00
+ ,1136029.00
+ ,10428.05
+ ,1241822.00
+ ,1498658.00
+ ,1623386.00
+ ,1741365.00)
+ ,dim=c(6
+ ,56)
+ ,dimnames=list(c('Y'
+ ,'DJIA'
+ ,'Y1'
+ ,'Y2'
+ ,'Y3'
+ ,'Y4')
+ ,1:56))
> y <- array(NA,dim=c(6,56),dimnames=list(c('Y','DJIA','Y1','Y2','Y3','Y4'),1:56))
> 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
Y DJIA Y1 Y2 Y3 Y4 M1 M2 M3 M4 M5 M6 M7 M8 M9
1 1469798 10467.48 1368839 1207763 1008380 989236 1 0 0 0 0 0 0 0 0
2 1498721 10274.97 1469798 1368839 1207763 1008380 0 1 0 0 0 0 0 0 0
3 1761769 10640.91 1498721 1469798 1368839 1207763 0 0 1 0 0 0 0 0 0
4 1653214 10481.60 1761769 1498721 1469798 1368839 0 0 0 1 0 0 0 0 0
5 1599104 10568.70 1653214 1761769 1498721 1469798 0 0 0 0 1 0 0 0 0
6 1421179 10440.07 1599104 1653214 1761769 1498721 0 0 0 0 0 1 0 0 0
7 1163995 10805.87 1421179 1599104 1653214 1761769 0 0 0 0 0 0 1 0 0
8 1037735 10717.50 1163995 1421179 1599104 1653214 0 0 0 0 0 0 0 1 0
9 1015407 10864.86 1037735 1163995 1421179 1599104 0 0 0 0 0 0 0 0 1
10 1039210 10993.41 1015407 1037735 1163995 1421179 0 0 0 0 0 0 0 0 0
11 1258049 11109.32 1039210 1015407 1037735 1163995 0 0 0 0 0 0 0 0 0
12 1469445 11367.14 1258049 1039210 1015407 1037735 0 0 0 0 0 0 0 0 0
13 1552346 11168.31 1469445 1258049 1039210 1015407 1 0 0 0 0 0 0 0 0
14 1549144 11150.22 1552346 1469445 1258049 1039210 0 1 0 0 0 0 0 0 0
15 1785895 11185.68 1549144 1552346 1469445 1258049 0 0 1 0 0 0 0 0 0
16 1662335 11381.15 1785895 1549144 1552346 1469445 0 0 0 1 0 0 0 0 0
17 1629440 11679.07 1662335 1785895 1549144 1552346 0 0 0 0 1 0 0 0 0
18 1467430 12080.73 1629440 1662335 1785895 1549144 0 0 0 0 0 1 0 0 0
19 1202209 12221.93 1467430 1629440 1662335 1785895 0 0 0 0 0 0 1 0 0
20 1076982 12463.15 1202209 1467430 1629440 1662335 0 0 0 0 0 0 0 1 0
21 1039367 12621.69 1076982 1202209 1467430 1629440 0 0 0 0 0 0 0 0 1
22 1063449 12268.63 1039367 1076982 1202209 1467430 0 0 0 0 0 0 0 0 0
23 1335135 12354.35 1063449 1039367 1076982 1202209 0 0 0 0 0 0 0 0 0
24 1491602 13062.91 1335135 1063449 1039367 1076982 0 0 0 0 0 0 0 0 0
25 1591972 13627.64 1491602 1335135 1063449 1039367 1 0 0 0 0 0 0 0 0
26 1641248 13408.62 1591972 1491602 1335135 1063449 0 1 0 0 0 0 0 0 0
27 1898849 13211.99 1641248 1591972 1491602 1335135 0 0 1 0 0 0 0 0 0
28 1798580 13357.74 1898849 1641248 1591972 1491602 0 0 0 1 0 0 0 0 0
29 1762444 13895.63 1798580 1898849 1641248 1591972 0 0 0 0 1 0 0 0 0
30 1622044 13930.01 1762444 1798580 1898849 1641248 0 0 0 0 0 1 0 0 0
31 1368955 13371.72 1622044 1762444 1798580 1898849 0 0 0 0 0 0 1 0 0
32 1262973 13264.82 1368955 1622044 1762444 1798580 0 0 0 0 0 0 0 1 0
33 1195650 12650.36 1262973 1368955 1622044 1762444 0 0 0 0 0 0 0 0 1
34 1269530 12266.39 1195650 1262973 1368955 1622044 0 0 0 0 0 0 0 0 0
35 1479279 12262.89 1269530 1195650 1262973 1368955 0 0 0 0 0 0 0 0 0
36 1607819 12820.13 1479279 1269530 1195650 1262973 0 0 0 0 0 0 0 0 0
37 1712466 12638.32 1607819 1479279 1269530 1195650 1 0 0 0 0 0 0 0 0
38 1721766 11350.01 1712466 1607819 1479279 1269530 0 1 0 0 0 0 0 0 0
39 1949843 11378.02 1721766 1712466 1607819 1479279 0 0 1 0 0 0 0 0 0
40 1821326 11543.55 1949843 1721766 1712466 1607819 0 0 0 1 0 0 0 0 0
41 1757802 10850.66 1821326 1949843 1721766 1712466 0 0 0 0 1 0 0 0 0
42 1590367 9325.01 1757802 1821326 1949843 1721766 0 0 0 0 0 1 0 0 0
43 1260647 8829.04 1590367 1757802 1821326 1949843 0 0 0 0 0 0 1 0 0
44 1149235 8776.39 1260647 1590367 1757802 1821326 0 0 0 0 0 0 0 1 0
45 1016367 8000.86 1149235 1260647 1590367 1757802 0 0 0 0 0 0 0 0 1
46 1027885 7062.93 1016367 1149235 1260647 1590367 0 0 0 0 0 0 0 0 0
47 1262159 7608.92 1027885 1016367 1149235 1260647 0 0 0 0 0 0 0 0 0
48 1520854 8168.12 1262159 1027885 1016367 1149235 0 0 0 0 0 0 0 0 0
49 1544144 8500.33 1520854 1262159 1027885 1016367 1 0 0 0 0 0 0 0 0
50 1564709 8447.00 1544144 1520854 1262159 1027885 0 1 0 0 0 0 0 0 0
51 1821776 9171.61 1564709 1544144 1520854 1262159 0 0 1 0 0 0 0 0 0
52 1741365 9496.28 1821776 1564709 1544144 1520854 0 0 0 1 0 0 0 0 0
53 1623386 9712.28 1741365 1821776 1564709 1544144 0 0 0 0 1 0 0 0 0
54 1498658 9712.73 1623386 1741365 1821776 1564709 0 0 0 0 0 1 0 0 0
55 1241822 10344.84 1498658 1623386 1741365 1821776 0 0 0 0 0 0 1 0 0
56 1136029 10428.05 1241822 1498658 1623386 1741365 0 0 0 0 0 0 0 1 0
M10 M11 t
1 0 0 1
2 0 0 2
3 0 0 3
4 0 0 4
5 0 0 5
6 0 0 6
7 0 0 7
8 0 0 8
9 0 0 9
10 1 0 10
11 0 1 11
12 0 0 12
13 0 0 13
14 0 0 14
15 0 0 15
16 0 0 16
17 0 0 17
18 0 0 18
19 0 0 19
20 0 0 20
21 0 0 21
22 1 0 22
23 0 1 23
24 0 0 24
25 0 0 25
26 0 0 26
27 0 0 27
28 0 0 28
29 0 0 29
30 0 0 30
31 0 0 31
32 0 0 32
33 0 0 33
34 1 0 34
35 0 1 35
36 0 0 36
37 0 0 37
38 0 0 38
39 0 0 39
40 0 0 40
41 0 0 41
42 0 0 42
43 0 0 43
44 0 0 44
45 0 0 45
46 1 0 46
47 0 1 47
48 0 0 48
49 0 0 49
50 0 0 50
51 0 0 51
52 0 0 52
53 0 0 53
54 0 0 54
55 0 0 55
56 0 0 56
> k <- length(x[1,])
> df <- as.data.frame(x)
> (mylm <- lm(df))
Call:
lm(formula = df)
Coefficients:
(Intercept) DJIA Y1 Y2 Y3 Y4
4.052e+05 1.429e+01 4.334e-01 1.173e-01 1.463e-03 1.982e-01
M1 M2 M3 M4 M5 M6
-2.055e+04 -5.884e+04 1.222e+05 -1.353e+05 -1.968e+05 -3.140e+05
M7 M8 M9 M10 M11 t
-5.624e+05 -5.206e+05 -4.897e+05 -3.771e+05 -9.877e+04 7.376e+02
> (mysum <- summary(mylm))
Call:
lm(formula = df)
Residuals:
Min 1Q Median 3Q Max
-49168 -13632 3103 12944 66725
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 4.052e+05 7.447e+04 5.441 3.32e-06 ***
DJIA 1.429e+01 3.482e+00 4.103 0.000208 ***
Y1 4.334e-01 1.532e-01 2.828 0.007438 **
Y2 1.173e-01 1.685e-01 0.696 0.490445
Y3 1.463e-03 1.667e-01 0.009 0.993043
Y4 1.982e-01 1.484e-01 1.336 0.189578
M1 -2.055e+04 4.176e+04 -0.492 0.625486
M2 -5.884e+04 6.151e+04 -0.957 0.344819
M3 1.222e+05 5.322e+04 2.297 0.027252 *
M4 -1.353e+05 4.655e+04 -2.907 0.006067 **
M5 -1.968e+05 6.301e+04 -3.124 0.003409 **
M6 -3.140e+05 6.773e+04 -4.636 4.11e-05 ***
M7 -5.624e+05 6.641e+04 -8.469 2.80e-10 ***
M8 -5.206e+05 8.579e+04 -6.068 4.59e-07 ***
M9 -4.897e+05 8.600e+04 -5.694 1.50e-06 ***
M10 -3.771e+05 7.266e+04 -5.190 7.32e-06 ***
M11 -9.877e+04 4.449e+04 -2.220 0.032469 *
t 7.376e+02 3.417e+02 2.158 0.037278 *
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 26130 on 38 degrees of freedom
Multiple R-squared: 0.9932, Adjusted R-squared: 0.9902
F-statistic: 326.8 on 17 and 38 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.023698668 0.04739734 0.9763013
[2,] 0.005923364 0.01184673 0.9940766
[3,] 0.021050890 0.04210178 0.9789491
[4,] 0.013074288 0.02614858 0.9869257
[5,] 0.021928555 0.04385711 0.9780714
[6,] 0.056588633 0.11317727 0.9434114
[7,] 0.069681170 0.13936234 0.9303188
[8,] 0.053361705 0.10672341 0.9466383
[9,] 0.026311678 0.05262336 0.9736883
[10,] 0.035711244 0.07142249 0.9642888
[11,] 0.024035299 0.04807060 0.9759647
[12,] 0.112290893 0.22458179 0.8877091
[13,] 0.236166917 0.47233383 0.7638331
[14,] 0.296619422 0.59323884 0.7033806
[15,] 0.498793084 0.99758617 0.5012069
> postscript(file="/var/www/html/rcomp/tmp/1y0h61292511666.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/rcomp/tmp/2y0h61292511666.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/rcomp/tmp/3y0h61292511666.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/rcomp/tmp/48ayr1292511666.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/rcomp/tmp/58ayr1292511666.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 = 56
Frequency = 1
1 2 3 4 5 6
2421.7296 4920.2495 16811.7004 17847.0039 19433.8779 -10151.3139
7 8 9 10 11 12
6561.5098 -7084.9944 32720.8535 1535.7156 -16921.4238 18714.6380
13 14 15 16 17 18
11375.6324 -19771.7526 -17345.0243 -31165.7016 1837.4501 -20437.6081
19 20 21 22 23 24
-12704.5832 -25468.2199 -4841.9573 -25505.8250 12579.0656 -36271.0907
25 26 27 28 29 30
-16415.9482 6525.7753 -2066.4344 3783.6820 14038.8806 6862.4942
31 32 33 34 35 36
23570.2811 22621.5232 15431.5579 51323.5510 8217.6232 -49168.3763
37 38 39 40 41 42
10813.1984 699.9605 -11486.9068 -11154.6468 4213.3031 15439.8400
43 44 45 46 47 48
-24552.1950 10313.2171 -43310.4541 -27353.4416 -3875.2650 66724.8290
49 50 51 52 53 54
-8194.6123 7625.7673 14086.6651 20689.6624 -39523.5118 8286.5878
55 56
7124.9873 -381.5260
> postscript(file="/var/www/html/rcomp/tmp/68ayr1292511666.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 = 56
Frequency = 1
lag(myerror, k = 1) myerror
0 2421.7296 NA
1 4920.2495 2421.7296
2 16811.7004 4920.2495
3 17847.0039 16811.7004
4 19433.8779 17847.0039
5 -10151.3139 19433.8779
6 6561.5098 -10151.3139
7 -7084.9944 6561.5098
8 32720.8535 -7084.9944
9 1535.7156 32720.8535
10 -16921.4238 1535.7156
11 18714.6380 -16921.4238
12 11375.6324 18714.6380
13 -19771.7526 11375.6324
14 -17345.0243 -19771.7526
15 -31165.7016 -17345.0243
16 1837.4501 -31165.7016
17 -20437.6081 1837.4501
18 -12704.5832 -20437.6081
19 -25468.2199 -12704.5832
20 -4841.9573 -25468.2199
21 -25505.8250 -4841.9573
22 12579.0656 -25505.8250
23 -36271.0907 12579.0656
24 -16415.9482 -36271.0907
25 6525.7753 -16415.9482
26 -2066.4344 6525.7753
27 3783.6820 -2066.4344
28 14038.8806 3783.6820
29 6862.4942 14038.8806
30 23570.2811 6862.4942
31 22621.5232 23570.2811
32 15431.5579 22621.5232
33 51323.5510 15431.5579
34 8217.6232 51323.5510
35 -49168.3763 8217.6232
36 10813.1984 -49168.3763
37 699.9605 10813.1984
38 -11486.9068 699.9605
39 -11154.6468 -11486.9068
40 4213.3031 -11154.6468
41 15439.8400 4213.3031
42 -24552.1950 15439.8400
43 10313.2171 -24552.1950
44 -43310.4541 10313.2171
45 -27353.4416 -43310.4541
46 -3875.2650 -27353.4416
47 66724.8290 -3875.2650
48 -8194.6123 66724.8290
49 7625.7673 -8194.6123
50 14086.6651 7625.7673
51 20689.6624 14086.6651
52 -39523.5118 20689.6624
53 8286.5878 -39523.5118
54 7124.9873 8286.5878
55 -381.5260 7124.9873
56 NA -381.5260
> dum1 <- dum[2:length(myerror),]
> dum1
lag(myerror, k = 1) myerror
[1,] 4920.2495 2421.7296
[2,] 16811.7004 4920.2495
[3,] 17847.0039 16811.7004
[4,] 19433.8779 17847.0039
[5,] -10151.3139 19433.8779
[6,] 6561.5098 -10151.3139
[7,] -7084.9944 6561.5098
[8,] 32720.8535 -7084.9944
[9,] 1535.7156 32720.8535
[10,] -16921.4238 1535.7156
[11,] 18714.6380 -16921.4238
[12,] 11375.6324 18714.6380
[13,] -19771.7526 11375.6324
[14,] -17345.0243 -19771.7526
[15,] -31165.7016 -17345.0243
[16,] 1837.4501 -31165.7016
[17,] -20437.6081 1837.4501
[18,] -12704.5832 -20437.6081
[19,] -25468.2199 -12704.5832
[20,] -4841.9573 -25468.2199
[21,] -25505.8250 -4841.9573
[22,] 12579.0656 -25505.8250
[23,] -36271.0907 12579.0656
[24,] -16415.9482 -36271.0907
[25,] 6525.7753 -16415.9482
[26,] -2066.4344 6525.7753
[27,] 3783.6820 -2066.4344
[28,] 14038.8806 3783.6820
[29,] 6862.4942 14038.8806
[30,] 23570.2811 6862.4942
[31,] 22621.5232 23570.2811
[32,] 15431.5579 22621.5232
[33,] 51323.5510 15431.5579
[34,] 8217.6232 51323.5510
[35,] -49168.3763 8217.6232
[36,] 10813.1984 -49168.3763
[37,] 699.9605 10813.1984
[38,] -11486.9068 699.9605
[39,] -11154.6468 -11486.9068
[40,] 4213.3031 -11154.6468
[41,] 15439.8400 4213.3031
[42,] -24552.1950 15439.8400
[43,] 10313.2171 -24552.1950
[44,] -43310.4541 10313.2171
[45,] -27353.4416 -43310.4541
[46,] -3875.2650 -27353.4416
[47,] 66724.8290 -3875.2650
[48,] -8194.6123 66724.8290
[49,] 7625.7673 -8194.6123
[50,] 14086.6651 7625.7673
[51,] 20689.6624 14086.6651
[52,] -39523.5118 20689.6624
[53,] 8286.5878 -39523.5118
[54,] 7124.9873 8286.5878
[55,] -381.5260 7124.9873
> z <- as.data.frame(dum1)
> z
lag(myerror, k = 1) myerror
1 4920.2495 2421.7296
2 16811.7004 4920.2495
3 17847.0039 16811.7004
4 19433.8779 17847.0039
5 -10151.3139 19433.8779
6 6561.5098 -10151.3139
7 -7084.9944 6561.5098
8 32720.8535 -7084.9944
9 1535.7156 32720.8535
10 -16921.4238 1535.7156
11 18714.6380 -16921.4238
12 11375.6324 18714.6380
13 -19771.7526 11375.6324
14 -17345.0243 -19771.7526
15 -31165.7016 -17345.0243
16 1837.4501 -31165.7016
17 -20437.6081 1837.4501
18 -12704.5832 -20437.6081
19 -25468.2199 -12704.5832
20 -4841.9573 -25468.2199
21 -25505.8250 -4841.9573
22 12579.0656 -25505.8250
23 -36271.0907 12579.0656
24 -16415.9482 -36271.0907
25 6525.7753 -16415.9482
26 -2066.4344 6525.7753
27 3783.6820 -2066.4344
28 14038.8806 3783.6820
29 6862.4942 14038.8806
30 23570.2811 6862.4942
31 22621.5232 23570.2811
32 15431.5579 22621.5232
33 51323.5510 15431.5579
34 8217.6232 51323.5510
35 -49168.3763 8217.6232
36 10813.1984 -49168.3763
37 699.9605 10813.1984
38 -11486.9068 699.9605
39 -11154.6468 -11486.9068
40 4213.3031 -11154.6468
41 15439.8400 4213.3031
42 -24552.1950 15439.8400
43 10313.2171 -24552.1950
44 -43310.4541 10313.2171
45 -27353.4416 -43310.4541
46 -3875.2650 -27353.4416
47 66724.8290 -3875.2650
48 -8194.6123 66724.8290
49 7625.7673 -8194.6123
50 14086.6651 7625.7673
51 20689.6624 14086.6651
52 -39523.5118 20689.6624
53 8286.5878 -39523.5118
54 7124.9873 8286.5878
55 -381.5260 7124.9873
> 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/711fc1292511666.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/rcomp/tmp/8usfx1292511666.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/rcomp/tmp/9usfx1292511666.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/rcomp/tmp/10mkei1292511666.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/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/11qkd61292511666.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/12bltc1292511666.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/137v931292511666.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/14tdp91292511666.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/15ew6w1292511666.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/16iem21292511666.tab")
+ }
>
> try(system("convert tmp/1y0h61292511666.ps tmp/1y0h61292511666.png",intern=TRUE))
character(0)
> try(system("convert tmp/2y0h61292511666.ps tmp/2y0h61292511666.png",intern=TRUE))
character(0)
> try(system("convert tmp/3y0h61292511666.ps tmp/3y0h61292511666.png",intern=TRUE))
character(0)
> try(system("convert tmp/48ayr1292511666.ps tmp/48ayr1292511666.png",intern=TRUE))
character(0)
> try(system("convert tmp/58ayr1292511666.ps tmp/58ayr1292511666.png",intern=TRUE))
character(0)
> try(system("convert tmp/68ayr1292511666.ps tmp/68ayr1292511666.png",intern=TRUE))
character(0)
> try(system("convert tmp/711fc1292511666.ps tmp/711fc1292511666.png",intern=TRUE))
character(0)
> try(system("convert tmp/8usfx1292511666.ps tmp/8usfx1292511666.png",intern=TRUE))
character(0)
> try(system("convert tmp/9usfx1292511666.ps tmp/9usfx1292511666.png",intern=TRUE))
character(0)
> try(system("convert tmp/10mkei1292511666.ps tmp/10mkei1292511666.png",intern=TRUE))
character(0)
>
>
> proc.time()
user system elapsed
2.395 1.699 8.430