R version 2.12.0 (2010-10-15)
Copyright (C) 2010 The R Foundation for Statistical Computing
ISBN 3-900051-07-0
Platform: i486-pc-linux-gnu (32-bit)
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(1368839.00
+ ,10192.51
+ ,1207763.00
+ ,1008380.00
+ ,989236.00
+ ,1469798.00
+ ,10467.48
+ ,1368839.00
+ ,1207763.00
+ ,1008380.00
+ ,1498721.00
+ ,10274.97
+ ,1469798.00
+ ,1368839.00
+ ,1207763.00
+ ,1761769.00
+ ,10640.91
+ ,1498721.00
+ ,1469798.00
+ ,1368839.00
+ ,1653214.00
+ ,10481.60
+ ,1761769.00
+ ,1498721.00
+ ,1469798.00
+ ,1599104.00
+ ,10568.70
+ ,1653214.00
+ ,1761769.00
+ ,1498721.00
+ ,1421179.00
+ ,10440.07
+ ,1599104.00
+ ,1653214.00
+ ,1761769.00
+ ,1163995.00
+ ,10805.87
+ ,1421179.00
+ ,1599104.00
+ ,1653214.00
+ ,1037735.00
+ ,10717.50
+ ,1163995.00
+ ,1421179.00
+ ,1599104.00
+ ,1015407.00
+ ,10864.86
+ ,1037735.00
+ ,1163995.00
+ ,1421179.00
+ ,1039210.00
+ ,10993.41
+ ,1015407.00
+ ,1037735.00
+ ,1163995.00
+ ,1258049.00
+ ,11109.32
+ ,1039210.00
+ ,1015407.00
+ ,1037735.00
+ ,1469445.00
+ ,11367.14
+ ,1258049.00
+ ,1039210.00
+ ,1015407.00
+ ,1552346.00
+ ,11168.31
+ ,1469445.00
+ ,1258049.00
+ ,1039210.00
+ ,1549144.00
+ ,11150.22
+ ,1552346.00
+ ,1469445.00
+ ,1258049.00
+ ,1785895.00
+ ,11185.68
+ ,1549144.00
+ ,1552346.00
+ ,1469445.00
+ ,1662335.00
+ ,11381.15
+ ,1785895.00
+ ,1549144.00
+ ,1552346.00
+ ,1629440.00
+ ,11679.07
+ ,1662335.00
+ ,1785895.00
+ ,1549144.00
+ ,1467430.00
+ ,12080.73
+ ,1629440.00
+ ,1662335.00
+ ,1785895.00
+ ,1202209.00
+ ,12221.93
+ ,1467430.00
+ ,1629440.00
+ ,1662335.00
+ ,1076982.00
+ ,12463.15
+ ,1202209.00
+ ,1467430.00
+ ,1629440.00
+ ,1039367.00
+ ,12621.69
+ ,1076982.00
+ ,1202209.00
+ ,1467430.00
+ ,1063449.00
+ ,12268.63
+ ,1039367.00
+ ,1076982.00
+ ,1202209.00
+ ,1335135.00
+ ,12354.35
+ ,1063449.00
+ ,1039367.00
+ ,1076982.00
+ ,1491602.00
+ ,13062.91
+ ,1335135.00
+ ,1063449.00
+ ,1039367.00
+ ,1591972.00
+ ,13627.64
+ ,1491602.00
+ ,1335135.00
+ ,1063449.00
+ ,1641248.00
+ ,13408.62
+ ,1591972.00
+ ,1491602.00
+ ,1335135.00
+ ,1898849.00
+ ,13211.99
+ ,1641248.00
+ ,1591972.00
+ ,1491602.00
+ ,1798580.00
+ ,13357.74
+ ,1898849.00
+ ,1641248.00
+ ,1591972.00
+ ,1762444.00
+ ,13895.63
+ ,1798580.00
+ ,1898849.00
+ ,1641248.00
+ ,1622044.00
+ ,13930.01
+ ,1762444.00
+ ,1798580.00
+ ,1898849.00
+ ,1368955.00
+ ,13371.72
+ ,1622044.00
+ ,1762444.00
+ ,1798580.00
+ ,1262973.00
+ ,13264.82
+ ,1368955.00
+ ,1622044.00
+ ,1762444.00
+ ,1195650.00
+ ,12650.36
+ ,1262973.00
+ ,1368955.00
+ ,1622044.00
+ ,1269530.00
+ ,12266.39
+ ,1195650.00
+ ,1262973.00
+ ,1368955.00
+ ,1479279.00
+ ,12262.89
+ ,1269530.00
+ ,1195650.00
+ ,1262973.00
+ ,1607819.00
+ ,12820.13
+ ,1479279.00
+ ,1269530.00
+ ,1195650.00
+ ,1712466.00
+ ,12638.32
+ ,1607819.00
+ ,1479279.00
+ ,1269530.00
+ ,1721766.00
+ ,11350.01
+ ,1712466.00
+ ,1607819.00
+ ,1479279.00
+ ,1949843.00
+ ,11378.02
+ ,1721766.00
+ ,1712466.00
+ ,1607819.00
+ ,1821326.00
+ ,11543.55
+ ,1949843.00
+ ,1721766.00
+ ,1712466.00
+ ,1757802.00
+ ,10850.66
+ ,1821326.00
+ ,1949843.00
+ ,1721766.00
+ ,1590367.00
+ ,9325.01
+ ,1757802.00
+ ,1821326.00
+ ,1949843.00
+ ,1260647.00
+ ,8829.04
+ ,1590367.00
+ ,1757802.00
+ ,1821326.00
+ ,1149235.00
+ ,8776.39
+ ,1260647.00
+ ,1590367.00
+ ,1757802.00
+ ,1016367.00
+ ,8000.86
+ ,1149235.00
+ ,1260647.00
+ ,1590367.00
+ ,1027885.00
+ ,7062.93
+ ,1016367.00
+ ,1149235.00
+ ,1260647.00
+ ,1262159.00
+ ,7608.92
+ ,1027885.00
+ ,1016367.00
+ ,1149235.00
+ ,1520854.00
+ ,8168.12
+ ,1262159.00
+ ,1027885.00
+ ,1016367.00
+ ,1544144.00
+ ,8500.33
+ ,1520854.00
+ ,1262159.00
+ ,1027885.00
+ ,1564709.00
+ ,8447.00
+ ,1544144.00
+ ,1520854.00
+ ,1262159.00
+ ,1821776.00
+ ,9171.61
+ ,1564709.00
+ ,1544144.00
+ ,1520854.00
+ ,1741365.00
+ ,9496.28
+ ,1821776.00
+ ,1564709.00
+ ,1544144.00
+ ,1623386.00
+ ,9712.28
+ ,1741365.00
+ ,1821776.00
+ ,1564709.00
+ ,1498658.00
+ ,9712.73
+ ,1623386.00
+ ,1741365.00
+ ,1821776.00
+ ,1241822.00
+ ,10344.84
+ ,1498658.00
+ ,1623386.00
+ ,1741365.00
+ ,1136029.00
+ ,10428.05
+ ,1241822.00
+ ,1498658.00
+ ,1623386.00)
+ ,dim=c(5
+ ,57)
+ ,dimnames=list(c('Y'
+ ,'DJIA'
+ ,'Y1'
+ ,'Y2'
+ ,'Y3')
+ ,1:57))
> y <- array(NA,dim=c(5,57),dimnames=list(c('Y','DJIA','Y1','Y2','Y3'),1:57))
> 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
> 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 M1 M2 M3 M4 M5 M6 M7 M8 M9 M10 M11
1 1368839 10192.51 1207763 1008380 989236 1 0 0 0 0 0 0 0 0 0 0
2 1469798 10467.48 1368839 1207763 1008380 0 1 0 0 0 0 0 0 0 0 0
3 1498721 10274.97 1469798 1368839 1207763 0 0 1 0 0 0 0 0 0 0 0
4 1761769 10640.91 1498721 1469798 1368839 0 0 0 1 0 0 0 0 0 0 0
5 1653214 10481.60 1761769 1498721 1469798 0 0 0 0 1 0 0 0 0 0 0
6 1599104 10568.70 1653214 1761769 1498721 0 0 0 0 0 1 0 0 0 0 0
7 1421179 10440.07 1599104 1653214 1761769 0 0 0 0 0 0 1 0 0 0 0
8 1163995 10805.87 1421179 1599104 1653214 0 0 0 0 0 0 0 1 0 0 0
9 1037735 10717.50 1163995 1421179 1599104 0 0 0 0 0 0 0 0 1 0 0
10 1015407 10864.86 1037735 1163995 1421179 0 0 0 0 0 0 0 0 0 1 0
11 1039210 10993.41 1015407 1037735 1163995 0 0 0 0 0 0 0 0 0 0 1
12 1258049 11109.32 1039210 1015407 1037735 0 0 0 0 0 0 0 0 0 0 0
13 1469445 11367.14 1258049 1039210 1015407 1 0 0 0 0 0 0 0 0 0 0
14 1552346 11168.31 1469445 1258049 1039210 0 1 0 0 0 0 0 0 0 0 0
15 1549144 11150.22 1552346 1469445 1258049 0 0 1 0 0 0 0 0 0 0 0
16 1785895 11185.68 1549144 1552346 1469445 0 0 0 1 0 0 0 0 0 0 0
17 1662335 11381.15 1785895 1549144 1552346 0 0 0 0 1 0 0 0 0 0 0
18 1629440 11679.07 1662335 1785895 1549144 0 0 0 0 0 1 0 0 0 0 0
19 1467430 12080.73 1629440 1662335 1785895 0 0 0 0 0 0 1 0 0 0 0
20 1202209 12221.93 1467430 1629440 1662335 0 0 0 0 0 0 0 1 0 0 0
21 1076982 12463.15 1202209 1467430 1629440 0 0 0 0 0 0 0 0 1 0 0
22 1039367 12621.69 1076982 1202209 1467430 0 0 0 0 0 0 0 0 0 1 0
23 1063449 12268.63 1039367 1076982 1202209 0 0 0 0 0 0 0 0 0 0 1
24 1335135 12354.35 1063449 1039367 1076982 0 0 0 0 0 0 0 0 0 0 0
25 1491602 13062.91 1335135 1063449 1039367 1 0 0 0 0 0 0 0 0 0 0
26 1591972 13627.64 1491602 1335135 1063449 0 1 0 0 0 0 0 0 0 0 0
27 1641248 13408.62 1591972 1491602 1335135 0 0 1 0 0 0 0 0 0 0 0
28 1898849 13211.99 1641248 1591972 1491602 0 0 0 1 0 0 0 0 0 0 0
29 1798580 13357.74 1898849 1641248 1591972 0 0 0 0 1 0 0 0 0 0 0
30 1762444 13895.63 1798580 1898849 1641248 0 0 0 0 0 1 0 0 0 0 0
31 1622044 13930.01 1762444 1798580 1898849 0 0 0 0 0 0 1 0 0 0 0
32 1368955 13371.72 1622044 1762444 1798580 0 0 0 0 0 0 0 1 0 0 0
33 1262973 13264.82 1368955 1622044 1762444 0 0 0 0 0 0 0 0 1 0 0
34 1195650 12650.36 1262973 1368955 1622044 0 0 0 0 0 0 0 0 0 1 0
35 1269530 12266.39 1195650 1262973 1368955 0 0 0 0 0 0 0 0 0 0 1
36 1479279 12262.89 1269530 1195650 1262973 0 0 0 0 0 0 0 0 0 0 0
37 1607819 12820.13 1479279 1269530 1195650 1 0 0 0 0 0 0 0 0 0 0
38 1712466 12638.32 1607819 1479279 1269530 0 1 0 0 0 0 0 0 0 0 0
39 1721766 11350.01 1712466 1607819 1479279 0 0 1 0 0 0 0 0 0 0 0
40 1949843 11378.02 1721766 1712466 1607819 0 0 0 1 0 0 0 0 0 0 0
41 1821326 11543.55 1949843 1721766 1712466 0 0 0 0 1 0 0 0 0 0 0
42 1757802 10850.66 1821326 1949843 1721766 0 0 0 0 0 1 0 0 0 0 0
43 1590367 9325.01 1757802 1821326 1949843 0 0 0 0 0 0 1 0 0 0 0
44 1260647 8829.04 1590367 1757802 1821326 0 0 0 0 0 0 0 1 0 0 0
45 1149235 8776.39 1260647 1590367 1757802 0 0 0 0 0 0 0 0 1 0 0
46 1016367 8000.86 1149235 1260647 1590367 0 0 0 0 0 0 0 0 0 1 0
47 1027885 7062.93 1016367 1149235 1260647 0 0 0 0 0 0 0 0 0 0 1
48 1262159 7608.92 1027885 1016367 1149235 0 0 0 0 0 0 0 0 0 0 0
49 1520854 8168.12 1262159 1027885 1016367 1 0 0 0 0 0 0 0 0 0 0
50 1544144 8500.33 1520854 1262159 1027885 0 1 0 0 0 0 0 0 0 0 0
51 1564709 8447.00 1544144 1520854 1262159 0 0 1 0 0 0 0 0 0 0 0
52 1821776 9171.61 1564709 1544144 1520854 0 0 0 1 0 0 0 0 0 0 0
53 1741365 9496.28 1821776 1564709 1544144 0 0 0 0 1 0 0 0 0 0 0
54 1623386 9712.28 1741365 1821776 1564709 0 0 0 0 0 1 0 0 0 0 0
55 1498658 9712.73 1623386 1741365 1821776 0 0 0 0 0 0 1 0 0 0 0
56 1241822 10344.84 1498658 1623386 1741365 0 0 0 0 0 0 0 1 0 0 0
57 1136029 10428.05 1241822 1498658 1623386 0 0 0 0 0 0 0 0 1 0 0
t
1 1
2 2
3 3
4 4
5 5
6 6
7 7
8 8
9 9
10 10
11 11
12 12
13 13
14 14
15 15
16 16
17 17
18 18
19 19
20 20
21 21
22 22
23 23
24 24
25 25
26 26
27 27
28 28
29 29
30 30
31 31
32 32
33 33
34 34
35 35
36 36
37 37
38 38
39 39
40 40
41 41
42 42
43 43
44 44
45 45
46 46
47 47
48 48
49 49
50 50
51 51
52 52
53 53
54 54
55 55
56 56
57 57
> k <- length(x[1,])
> df <- as.data.frame(x)
> (mylm <- lm(df))
Call:
lm(formula = df)
Coefficients:
(Intercept) DJIA Y1 Y2 Y3 M1
3.460e+05 1.287e+01 4.882e-01 1.490e-01 1.132e-01 6.351e+04
M2 M3 M4 M5 M6 M7
1.642e+04 -5.207e+04 1.499e+05 -9.466e+04 -1.441e+05 -2.785e+05
M8 M9 M10 M11 t
-4.552e+05 -4.083e+05 -3.551e+05 -2.371e+05 8.095e+02
> (mysum <- summary(mylm))
Call:
lm(formula = df)
Residuals:
Min 1Q Median 3Q Max
-45751 -14896 2488 10877 82291
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 3.460e+05 5.924e+04 5.840 7.94e-07 ***
DJIA 1.287e+01 3.242e+00 3.969 0.000292 ***
Y1 4.882e-01 1.518e-01 3.216 0.002572 **
Y2 1.490e-01 1.670e-01 0.892 0.377758
Y3 1.132e-01 1.399e-01 0.809 0.423240
M1 6.351e+04 3.990e+04 1.592 0.119335
M2 1.642e+04 6.005e+04 0.273 0.785961
M3 -5.207e+04 5.385e+04 -0.967 0.339379
M4 1.499e+05 4.465e+04 3.358 0.001735 **
M5 -9.466e+04 6.402e+04 -1.479 0.147061
M6 -1.441e+05 6.536e+04 -2.205 0.033281 *
M7 -2.785e+05 4.729e+04 -5.889 6.77e-07 ***
M8 -4.552e+05 4.381e+04 -10.391 6.31e-13 ***
M9 -4.083e+05 5.209e+04 -7.837 1.32e-09 ***
M10 -3.551e+05 4.504e+04 -7.885 1.13e-09 ***
M11 -2.371e+05 2.527e+04 -9.383 1.17e-11 ***
t 8.095e+02 3.380e+02 2.395 0.021403 *
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 26450 on 40 degrees of freedom
Multiple R-squared: 0.9927, Adjusted R-squared: 0.9898
F-statistic: 339.3 on 16 and 40 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.23184094 0.46368188 0.7681591
[2,] 0.11714990 0.23429980 0.8828501
[3,] 0.05163659 0.10327317 0.9483634
[4,] 0.02601421 0.05202842 0.9739858
[5,] 0.02912104 0.05824208 0.9708790
[6,] 0.02662540 0.05325080 0.9733746
[7,] 0.02130133 0.04260266 0.9786987
[8,] 0.06368902 0.12737804 0.9363110
[9,] 0.04012978 0.08025956 0.9598702
[10,] 0.03334100 0.06668201 0.9666590
[11,] 0.01633648 0.03267297 0.9836635
[12,] 0.02271753 0.04543506 0.9772825
[13,] 0.01317316 0.02634631 0.9868268
[14,] 0.08327631 0.16655263 0.9167237
[15,] 0.15861108 0.31722215 0.8413889
[16,] 0.39596238 0.79192477 0.6040376
[17,] 0.31648992 0.63297983 0.6835101
[18,] 0.33560918 0.67121836 0.6643908
> postscript(file="/var/www/rcomp/tmp/1wvm31292511531.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/rcomp/tmp/2wvm31292511531.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/rcomp/tmp/3p53o1292511531.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/rcomp/tmp/4p53o1292511531.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/rcomp/tmp/5p53o1292511531.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 = 57
Frequency = 1
1 2 3 4 5 6
-24383.7528 8823.4549 12052.4961 20214.1594 13321.8352 17244.6136
7 8 9 10 11 12
-12590.2157 8600.4509 -6105.6538 35784.0432 -2026.2469 -16614.8044
13 14 15 16 17 18
19292.8247 12547.7154 -19475.8021 -20682.4675 -27477.2905 9826.8330
19 20 21 22 23 24
-16060.0819 -9247.9747 -28010.4706 -2659.9708 -25776.5625 14894.5876
25 26 27 28 29 30
-34035.1215 -14225.3949 2487.9421 3114.4413 276.8105 10836.5850
31 32 33 34 35 36
7037.3734 22278.3945 18464.8471 10404.5890 47745.7413 5565.7398
37 38 39 40 41 42
-43164.8907 7750.2252 7328.1641 -2420.8778 -13877.1008 7842.6375
43 44 45 46 47 48
17999.5167 -23717.2375 10877.2447 -43528.6614 -19942.9320 -3845.5230
49 50 51 52 53 54
82290.9402 -14896.0006 -2392.8002 -225.2555 27755.7454 -45750.6691
55 56 57
3613.4075 2086.3668 4774.0326
> postscript(file="/var/www/rcomp/tmp/6iwlr1292511531.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 = 57
Frequency = 1
lag(myerror, k = 1) myerror
0 -24383.7528 NA
1 8823.4549 -24383.7528
2 12052.4961 8823.4549
3 20214.1594 12052.4961
4 13321.8352 20214.1594
5 17244.6136 13321.8352
6 -12590.2157 17244.6136
7 8600.4509 -12590.2157
8 -6105.6538 8600.4509
9 35784.0432 -6105.6538
10 -2026.2469 35784.0432
11 -16614.8044 -2026.2469
12 19292.8247 -16614.8044
13 12547.7154 19292.8247
14 -19475.8021 12547.7154
15 -20682.4675 -19475.8021
16 -27477.2905 -20682.4675
17 9826.8330 -27477.2905
18 -16060.0819 9826.8330
19 -9247.9747 -16060.0819
20 -28010.4706 -9247.9747
21 -2659.9708 -28010.4706
22 -25776.5625 -2659.9708
23 14894.5876 -25776.5625
24 -34035.1215 14894.5876
25 -14225.3949 -34035.1215
26 2487.9421 -14225.3949
27 3114.4413 2487.9421
28 276.8105 3114.4413
29 10836.5850 276.8105
30 7037.3734 10836.5850
31 22278.3945 7037.3734
32 18464.8471 22278.3945
33 10404.5890 18464.8471
34 47745.7413 10404.5890
35 5565.7398 47745.7413
36 -43164.8907 5565.7398
37 7750.2252 -43164.8907
38 7328.1641 7750.2252
39 -2420.8778 7328.1641
40 -13877.1008 -2420.8778
41 7842.6375 -13877.1008
42 17999.5167 7842.6375
43 -23717.2375 17999.5167
44 10877.2447 -23717.2375
45 -43528.6614 10877.2447
46 -19942.9320 -43528.6614
47 -3845.5230 -19942.9320
48 82290.9402 -3845.5230
49 -14896.0006 82290.9402
50 -2392.8002 -14896.0006
51 -225.2555 -2392.8002
52 27755.7454 -225.2555
53 -45750.6691 27755.7454
54 3613.4075 -45750.6691
55 2086.3668 3613.4075
56 4774.0326 2086.3668
57 NA 4774.0326
> dum1 <- dum[2:length(myerror),]
> dum1
lag(myerror, k = 1) myerror
[1,] 8823.4549 -24383.7528
[2,] 12052.4961 8823.4549
[3,] 20214.1594 12052.4961
[4,] 13321.8352 20214.1594
[5,] 17244.6136 13321.8352
[6,] -12590.2157 17244.6136
[7,] 8600.4509 -12590.2157
[8,] -6105.6538 8600.4509
[9,] 35784.0432 -6105.6538
[10,] -2026.2469 35784.0432
[11,] -16614.8044 -2026.2469
[12,] 19292.8247 -16614.8044
[13,] 12547.7154 19292.8247
[14,] -19475.8021 12547.7154
[15,] -20682.4675 -19475.8021
[16,] -27477.2905 -20682.4675
[17,] 9826.8330 -27477.2905
[18,] -16060.0819 9826.8330
[19,] -9247.9747 -16060.0819
[20,] -28010.4706 -9247.9747
[21,] -2659.9708 -28010.4706
[22,] -25776.5625 -2659.9708
[23,] 14894.5876 -25776.5625
[24,] -34035.1215 14894.5876
[25,] -14225.3949 -34035.1215
[26,] 2487.9421 -14225.3949
[27,] 3114.4413 2487.9421
[28,] 276.8105 3114.4413
[29,] 10836.5850 276.8105
[30,] 7037.3734 10836.5850
[31,] 22278.3945 7037.3734
[32,] 18464.8471 22278.3945
[33,] 10404.5890 18464.8471
[34,] 47745.7413 10404.5890
[35,] 5565.7398 47745.7413
[36,] -43164.8907 5565.7398
[37,] 7750.2252 -43164.8907
[38,] 7328.1641 7750.2252
[39,] -2420.8778 7328.1641
[40,] -13877.1008 -2420.8778
[41,] 7842.6375 -13877.1008
[42,] 17999.5167 7842.6375
[43,] -23717.2375 17999.5167
[44,] 10877.2447 -23717.2375
[45,] -43528.6614 10877.2447
[46,] -19942.9320 -43528.6614
[47,] -3845.5230 -19942.9320
[48,] 82290.9402 -3845.5230
[49,] -14896.0006 82290.9402
[50,] -2392.8002 -14896.0006
[51,] -225.2555 -2392.8002
[52,] 27755.7454 -225.2555
[53,] -45750.6691 27755.7454
[54,] 3613.4075 -45750.6691
[55,] 2086.3668 3613.4075
[56,] 4774.0326 2086.3668
> z <- as.data.frame(dum1)
> z
lag(myerror, k = 1) myerror
1 8823.4549 -24383.7528
2 12052.4961 8823.4549
3 20214.1594 12052.4961
4 13321.8352 20214.1594
5 17244.6136 13321.8352
6 -12590.2157 17244.6136
7 8600.4509 -12590.2157
8 -6105.6538 8600.4509
9 35784.0432 -6105.6538
10 -2026.2469 35784.0432
11 -16614.8044 -2026.2469
12 19292.8247 -16614.8044
13 12547.7154 19292.8247
14 -19475.8021 12547.7154
15 -20682.4675 -19475.8021
16 -27477.2905 -20682.4675
17 9826.8330 -27477.2905
18 -16060.0819 9826.8330
19 -9247.9747 -16060.0819
20 -28010.4706 -9247.9747
21 -2659.9708 -28010.4706
22 -25776.5625 -2659.9708
23 14894.5876 -25776.5625
24 -34035.1215 14894.5876
25 -14225.3949 -34035.1215
26 2487.9421 -14225.3949
27 3114.4413 2487.9421
28 276.8105 3114.4413
29 10836.5850 276.8105
30 7037.3734 10836.5850
31 22278.3945 7037.3734
32 18464.8471 22278.3945
33 10404.5890 18464.8471
34 47745.7413 10404.5890
35 5565.7398 47745.7413
36 -43164.8907 5565.7398
37 7750.2252 -43164.8907
38 7328.1641 7750.2252
39 -2420.8778 7328.1641
40 -13877.1008 -2420.8778
41 7842.6375 -13877.1008
42 17999.5167 7842.6375
43 -23717.2375 17999.5167
44 10877.2447 -23717.2375
45 -43528.6614 10877.2447
46 -19942.9320 -43528.6614
47 -3845.5230 -19942.9320
48 82290.9402 -3845.5230
49 -14896.0006 82290.9402
50 -2392.8002 -14896.0006
51 -225.2555 -2392.8002
52 27755.7454 -225.2555
53 -45750.6691 27755.7454
54 3613.4075 -45750.6691
55 2086.3668 3613.4075
56 4774.0326 2086.3668
> 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/rcomp/tmp/7sn2c1292511531.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/rcomp/tmp/8sn2c1292511531.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/rcomp/tmp/9sn2c1292511531.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/rcomp/tmp/103ejw1292511531.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/rcomp/createtable file can be downloaded at http://www.wessa.net/cretab
> load(file="/var/www/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/rcomp/tmp/11pf021292511531.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/rcomp/tmp/12sfg81292511531.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/rcomp/tmp/13bsl51292511531.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/rcomp/tmp/14sqv51292511531.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/rcomp/tmp/15v8bt1292511531.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/rcomp/tmp/16gr9z1292511531.tab")
+ }
>
> try(system("convert tmp/1wvm31292511531.ps tmp/1wvm31292511531.png",intern=TRUE))
character(0)
> try(system("convert tmp/2wvm31292511531.ps tmp/2wvm31292511531.png",intern=TRUE))
character(0)
> try(system("convert tmp/3p53o1292511531.ps tmp/3p53o1292511531.png",intern=TRUE))
character(0)
> try(system("convert tmp/4p53o1292511531.ps tmp/4p53o1292511531.png",intern=TRUE))
character(0)
> try(system("convert tmp/5p53o1292511531.ps tmp/5p53o1292511531.png",intern=TRUE))
character(0)
> try(system("convert tmp/6iwlr1292511531.ps tmp/6iwlr1292511531.png",intern=TRUE))
character(0)
> try(system("convert tmp/7sn2c1292511531.ps tmp/7sn2c1292511531.png",intern=TRUE))
character(0)
> try(system("convert tmp/8sn2c1292511531.ps tmp/8sn2c1292511531.png",intern=TRUE))
character(0)
> try(system("convert tmp/9sn2c1292511531.ps tmp/9sn2c1292511531.png",intern=TRUE))
character(0)
> try(system("convert tmp/103ejw1292511531.ps tmp/103ejw1292511531.png",intern=TRUE))
character(0)
>
>
> proc.time()
user system elapsed
3.080 1.660 4.725