R version 2.8.0 (2008-10-20)
Copyright (C) 2008 The R Foundation for Statistical Computing
ISBN 3-900051-07-0
R is free software and comes with ABSOLUTELY NO WARRANTY.
You are welcome to redistribute it under certain conditions.
Type 'license()' or 'licence()' for distribution details.
Natural language support but running in an English locale
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(1207763.00
+ ,10503.76
+ ,1008380.00
+ ,989236.00
+ ,1368839.00
+ ,10192.51
+ ,1207763.00
+ ,1008380.00
+ ,1469798.00
+ ,10467.48
+ ,1368839.00
+ ,1207763.00
+ ,1498721.00
+ ,10274.97
+ ,1469798.00
+ ,1368839.00
+ ,1761769.00
+ ,10640.91
+ ,1498721.00
+ ,1469798.00
+ ,1653214.00
+ ,10481.60
+ ,1761769.00
+ ,1498721.00
+ ,1599104.00
+ ,10568.70
+ ,1653214.00
+ ,1761769.00
+ ,1421179.00
+ ,10440.07
+ ,1599104.00
+ ,1653214.00
+ ,1163995.00
+ ,10805.87
+ ,1421179.00
+ ,1599104.00
+ ,1037735.00
+ ,10717.50
+ ,1163995.00
+ ,1421179.00
+ ,1015407.00
+ ,10864.86
+ ,1037735.00
+ ,1163995.00
+ ,1039210.00
+ ,10993.41
+ ,1015407.00
+ ,1037735.00
+ ,1258049.00
+ ,11109.32
+ ,1039210.00
+ ,1015407.00
+ ,1469445.00
+ ,11367.14
+ ,1258049.00
+ ,1039210.00
+ ,1552346.00
+ ,11168.31
+ ,1469445.00
+ ,1258049.00
+ ,1549144.00
+ ,11150.22
+ ,1552346.00
+ ,1469445.00
+ ,1785895.00
+ ,11185.68
+ ,1549144.00
+ ,1552346.00
+ ,1662335.00
+ ,11381.15
+ ,1785895.00
+ ,1549144.00
+ ,1629440.00
+ ,11679.07
+ ,1662335.00
+ ,1785895.00
+ ,1467430.00
+ ,12080.73
+ ,1629440.00
+ ,1662335.00
+ ,1202209.00
+ ,12221.93
+ ,1467430.00
+ ,1629440.00
+ ,1076982.00
+ ,12463.15
+ ,1202209.00
+ ,1467430.00
+ ,1039367.00
+ ,12621.69
+ ,1076982.00
+ ,1202209.00
+ ,1063449.00
+ ,12268.63
+ ,1039367.00
+ ,1076982.00
+ ,1335135.00
+ ,12354.35
+ ,1063449.00
+ ,1039367.00
+ ,1491602.00
+ ,13062.91
+ ,1335135.00
+ ,1063449.00
+ ,1591972.00
+ ,13627.64
+ ,1491602.00
+ ,1335135.00
+ ,1641248.00
+ ,13408.62
+ ,1591972.00
+ ,1491602.00
+ ,1898849.00
+ ,13211.99
+ ,1641248.00
+ ,1591972.00
+ ,1798580.00
+ ,13357.74
+ ,1898849.00
+ ,1641248.00
+ ,1762444.00
+ ,13895.63
+ ,1798580.00
+ ,1898849.00
+ ,1622044.00
+ ,13930.01
+ ,1762444.00
+ ,1798580.00
+ ,1368955.00
+ ,13371.72
+ ,1622044.00
+ ,1762444.00
+ ,1262973.00
+ ,13264.82
+ ,1368955.00
+ ,1622044.00
+ ,1195650.00
+ ,12650.36
+ ,1262973.00
+ ,1368955.00
+ ,1269530.00
+ ,12266.39
+ ,1195650.00
+ ,1262973.00
+ ,1479279.00
+ ,12262.89
+ ,1269530.00
+ ,1195650.00
+ ,1607819.00
+ ,12820.13
+ ,1479279.00
+ ,1269530.00
+ ,1712466.00
+ ,12638.32
+ ,1607819.00
+ ,1479279.00
+ ,1721766.00
+ ,11350.01
+ ,1712466.00
+ ,1607819.00
+ ,1949843.00
+ ,11378.02
+ ,1721766.00
+ ,1712466.00
+ ,1821326.00
+ ,11543.55
+ ,1949843.00
+ ,1721766.00
+ ,1757802.00
+ ,10850.66
+ ,1821326.00
+ ,1949843.00
+ ,1590367.00
+ ,9325.01
+ ,1757802.00
+ ,1821326.00
+ ,1260647.00
+ ,8829.04
+ ,1590367.00
+ ,1757802.00
+ ,1149235.00
+ ,8776.39
+ ,1260647.00
+ ,1590367.00
+ ,1016367.00
+ ,8000.86
+ ,1149235.00
+ ,1260647.00
+ ,1027885.00
+ ,7062.93
+ ,1016367.00
+ ,1149235.00
+ ,1262159.00
+ ,7608.92
+ ,1027885.00
+ ,1016367.00
+ ,1520854.00
+ ,8168.12
+ ,1262159.00
+ ,1027885.00
+ ,1544144.00
+ ,8500.33
+ ,1520854.00
+ ,1262159.00
+ ,1564709.00
+ ,8447.00
+ ,1544144.00
+ ,1520854.00
+ ,1821776.00
+ ,9171.61
+ ,1564709.00
+ ,1544144.00
+ ,1741365.00
+ ,9496.28
+ ,1821776.00
+ ,1564709.00
+ ,1623386.00
+ ,9712.28
+ ,1741365.00
+ ,1821776.00
+ ,1498658.00
+ ,9712.73
+ ,1623386.00
+ ,1741365.00
+ ,1241822.00
+ ,10344.84
+ ,1498658.00
+ ,1623386.00
+ ,1136029.00
+ ,10428.05
+ ,1241822.00
+ ,1498658.00)
+ ,dim=c(4
+ ,58)
+ ,dimnames=list(c('Y'
+ ,'DJIA'
+ ,'Y1'
+ ,'Y2')
+ ,1:58))
> y <- array(NA,dim=c(4,58),dimnames=list(c('Y','DJIA','Y1','Y2'),1:58))
> 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 M1 M2 M3 M4 M5 M6 M7 M8 M9 M10 M11 t
1 1207763 10503.76 1008380 989236 1 0 0 0 0 0 0 0 0 0 0 1
2 1368839 10192.51 1207763 1008380 0 1 0 0 0 0 0 0 0 0 0 2
3 1469798 10467.48 1368839 1207763 0 0 1 0 0 0 0 0 0 0 0 3
4 1498721 10274.97 1469798 1368839 0 0 0 1 0 0 0 0 0 0 0 4
5 1761769 10640.91 1498721 1469798 0 0 0 0 1 0 0 0 0 0 0 5
6 1653214 10481.60 1761769 1498721 0 0 0 0 0 1 0 0 0 0 0 6
7 1599104 10568.70 1653214 1761769 0 0 0 0 0 0 1 0 0 0 0 7
8 1421179 10440.07 1599104 1653214 0 0 0 0 0 0 0 1 0 0 0 8
9 1163995 10805.87 1421179 1599104 0 0 0 0 0 0 0 0 1 0 0 9
10 1037735 10717.50 1163995 1421179 0 0 0 0 0 0 0 0 0 1 0 10
11 1015407 10864.86 1037735 1163995 0 0 0 0 0 0 0 0 0 0 1 11
12 1039210 10993.41 1015407 1037735 0 0 0 0 0 0 0 0 0 0 0 12
13 1258049 11109.32 1039210 1015407 1 0 0 0 0 0 0 0 0 0 0 13
14 1469445 11367.14 1258049 1039210 0 1 0 0 0 0 0 0 0 0 0 14
15 1552346 11168.31 1469445 1258049 0 0 1 0 0 0 0 0 0 0 0 15
16 1549144 11150.22 1552346 1469445 0 0 0 1 0 0 0 0 0 0 0 16
17 1785895 11185.68 1549144 1552346 0 0 0 0 1 0 0 0 0 0 0 17
18 1662335 11381.15 1785895 1549144 0 0 0 0 0 1 0 0 0 0 0 18
19 1629440 11679.07 1662335 1785895 0 0 0 0 0 0 1 0 0 0 0 19
20 1467430 12080.73 1629440 1662335 0 0 0 0 0 0 0 1 0 0 0 20
21 1202209 12221.93 1467430 1629440 0 0 0 0 0 0 0 0 1 0 0 21
22 1076982 12463.15 1202209 1467430 0 0 0 0 0 0 0 0 0 1 0 22
23 1039367 12621.69 1076982 1202209 0 0 0 0 0 0 0 0 0 0 1 23
24 1063449 12268.63 1039367 1076982 0 0 0 0 0 0 0 0 0 0 0 24
25 1335135 12354.35 1063449 1039367 1 0 0 0 0 0 0 0 0 0 0 25
26 1491602 13062.91 1335135 1063449 0 1 0 0 0 0 0 0 0 0 0 26
27 1591972 13627.64 1491602 1335135 0 0 1 0 0 0 0 0 0 0 0 27
28 1641248 13408.62 1591972 1491602 0 0 0 1 0 0 0 0 0 0 0 28
29 1898849 13211.99 1641248 1591972 0 0 0 0 1 0 0 0 0 0 0 29
30 1798580 13357.74 1898849 1641248 0 0 0 0 0 1 0 0 0 0 0 30
31 1762444 13895.63 1798580 1898849 0 0 0 0 0 0 1 0 0 0 0 31
32 1622044 13930.01 1762444 1798580 0 0 0 0 0 0 0 1 0 0 0 32
33 1368955 13371.72 1622044 1762444 0 0 0 0 0 0 0 0 1 0 0 33
34 1262973 13264.82 1368955 1622044 0 0 0 0 0 0 0 0 0 1 0 34
35 1195650 12650.36 1262973 1368955 0 0 0 0 0 0 0 0 0 0 1 35
36 1269530 12266.39 1195650 1262973 0 0 0 0 0 0 0 0 0 0 0 36
37 1479279 12262.89 1269530 1195650 1 0 0 0 0 0 0 0 0 0 0 37
38 1607819 12820.13 1479279 1269530 0 1 0 0 0 0 0 0 0 0 0 38
39 1712466 12638.32 1607819 1479279 0 0 1 0 0 0 0 0 0 0 0 39
40 1721766 11350.01 1712466 1607819 0 0 0 1 0 0 0 0 0 0 0 40
41 1949843 11378.02 1721766 1712466 0 0 0 0 1 0 0 0 0 0 0 41
42 1821326 11543.55 1949843 1721766 0 0 0 0 0 1 0 0 0 0 0 42
43 1757802 10850.66 1821326 1949843 0 0 0 0 0 0 1 0 0 0 0 43
44 1590367 9325.01 1757802 1821326 0 0 0 0 0 0 0 1 0 0 0 44
45 1260647 8829.04 1590367 1757802 0 0 0 0 0 0 0 0 1 0 0 45
46 1149235 8776.39 1260647 1590367 0 0 0 0 0 0 0 0 0 1 0 46
47 1016367 8000.86 1149235 1260647 0 0 0 0 0 0 0 0 0 0 1 47
48 1027885 7062.93 1016367 1149235 0 0 0 0 0 0 0 0 0 0 0 48
49 1262159 7608.92 1027885 1016367 1 0 0 0 0 0 0 0 0 0 0 49
50 1520854 8168.12 1262159 1027885 0 1 0 0 0 0 0 0 0 0 0 50
51 1544144 8500.33 1520854 1262159 0 0 1 0 0 0 0 0 0 0 0 51
52 1564709 8447.00 1544144 1520854 0 0 0 1 0 0 0 0 0 0 0 52
53 1821776 9171.61 1564709 1544144 0 0 0 0 1 0 0 0 0 0 0 53
54 1741365 9496.28 1821776 1564709 0 0 0 0 0 1 0 0 0 0 0 54
55 1623386 9712.28 1741365 1821776 0 0 0 0 0 0 1 0 0 0 0 55
56 1498658 9712.73 1623386 1741365 0 0 0 0 0 0 0 1 0 0 0 56
57 1241822 10344.84 1498658 1623386 0 0 0 0 0 0 0 0 1 0 0 57
58 1136029 10428.05 1241822 1498658 0 0 0 0 0 0 0 0 0 1 0 58
> k <- length(x[1,])
> df <- as.data.frame(x)
> (mylm <- lm(df))
Call:
lm(formula = df)
Coefficients:
(Intercept) DJIA Y1 Y2 M1 M2
1.297e+05 1.244e+01 5.248e-01 2.226e-01 2.214e+05 2.735e+05
M3 M4 M5 M6 M7 M8
2.065e+05 1.469e+05 3.628e+05 1.169e+05 5.555e+04 -4.073e+04
M9 M10 M11 t
-2.196e+05 -1.582e+05 -9.747e+04 8.746e+02
> (mysum <- summary(mylm))
Call:
lm(formula = df)
Residuals:
Min 1Q Median 3Q Max
-50068 -16521 2498 12851 81172
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 1.297e+05 5.764e+04 2.249 0.029781 *
DJIA 1.244e+01 3.079e+00 4.039 0.000224 ***
Y1 5.248e-01 1.458e-01 3.599 0.000837 ***
Y2 2.226e-01 1.323e-01 1.682 0.099971 .
M1 2.214e+05 2.175e+04 10.177 6.62e-13 ***
M2 2.735e+05 4.460e+04 6.134 2.55e-07 ***
M3 2.065e+05 4.597e+04 4.492 5.43e-05 ***
M4 1.469e+05 4.338e+04 3.387 0.001543 **
M5 3.628e+05 4.161e+04 8.721 5.65e-11 ***
M6 1.169e+05 6.785e+04 1.724 0.092127 .
M7 5.555e+04 4.981e+04 1.115 0.271119
M8 -4.073e+04 4.526e+04 -0.900 0.373330
M9 -2.196e+05 3.792e+04 -5.791 7.93e-07 ***
M10 -1.582e+05 3.526e+04 -4.487 5.52e-05 ***
M11 -9.747e+04 2.049e+04 -4.757 2.33e-05 ***
t 8.746e+02 3.302e+02 2.649 0.011330 *
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Residual standard error: 26380 on 42 degrees of freedom
Multiple R-squared: 0.9925, Adjusted R-squared: 0.9898
F-statistic: 369.7 on 15 and 42 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.127475651 0.25495130 0.8725243
[2,] 0.109202513 0.21840503 0.8907975
[3,] 0.065791359 0.13158272 0.9342086
[4,] 0.029288467 0.05857693 0.9707115
[5,] 0.015325170 0.03065034 0.9846748
[6,] 0.006798846 0.01359769 0.9932012
[7,] 0.009819041 0.01963808 0.9901810
[8,] 0.015851229 0.03170246 0.9841488
[9,] 0.025054441 0.05010888 0.9749456
[10,] 0.031888859 0.06377772 0.9681111
[11,] 0.023212224 0.04642445 0.9767878
[12,] 0.026889029 0.05377806 0.9731110
[13,] 0.015715108 0.03143022 0.9842849
[14,] 0.029279989 0.05855998 0.9707200
[15,] 0.023422532 0.04684506 0.9765775
[16,] 0.164931241 0.32986248 0.8350688
[17,] 0.193152001 0.38630400 0.8068480
[18,] 0.194751375 0.38950275 0.8052486
[19,] 0.278883599 0.55776720 0.7211164
[20,] 0.480762535 0.96152507 0.5192375
[21,] 0.341184747 0.68236949 0.6588153
> postscript(file="/var/www/html/freestat/rcomp/tmp/1v8bt1292511296.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/freestat/rcomp/tmp/2v8bt1292511296.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/freestat/rcomp/tmp/3v8bt1292511296.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/freestat/rcomp/tmp/46hse1292511296.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/freestat/rcomp/tmp/56hse1292511296.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 = 58
Frequency = 1
1 2 3 4 5 6
-24171.8462 -21147.7309 13654.7690 14802.9445 18876.8096 12837.0027
7 8 9 10 11 12
16599.6766 -11768.5107 9896.1145 -2917.2423 34800.7368 -1524.4127
13 14 15 16 17 18
-13916.8610 21102.4566 12999.7488 -21868.1109 -19102.4290 -23608.1563
19 20 21 22 23 24
12475.4609 -14366.9931 -11020.4817 -26223.6199 -2685.0977 -24948.9035
25 26 27 28 29 30
19136.9475 -34175.6075 -17238.4220 5927.3334 999.5841 -2218.9730
31 32 33 34 35 36
10775.3140 6627.6100 20184.7555 17379.0862 8022.1898 47249.0959
37 38 39 40 41 42
10984.0783 -46950.6737 11988.9330 12447.0404 -4769.7976 -12090.3169
43 44 45 46 47 48
10217.9553 19096.4171 -24468.3063 12856.1676 -40137.8288 -20775.7797
49 50 51 52 53 54
7967.6814 81171.5554 -21405.0288 -11309.2074 3995.8328 25080.4435
55 56 57 58
-50068.4068 411.4767 5407.9179 -1094.3915
> postscript(file="/var/www/html/freestat/rcomp/tmp/66hse1292511296.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 = 58
Frequency = 1
lag(myerror, k = 1) myerror
0 -24171.8462 NA
1 -21147.7309 -24171.8462
2 13654.7690 -21147.7309
3 14802.9445 13654.7690
4 18876.8096 14802.9445
5 12837.0027 18876.8096
6 16599.6766 12837.0027
7 -11768.5107 16599.6766
8 9896.1145 -11768.5107
9 -2917.2423 9896.1145
10 34800.7368 -2917.2423
11 -1524.4127 34800.7368
12 -13916.8610 -1524.4127
13 21102.4566 -13916.8610
14 12999.7488 21102.4566
15 -21868.1109 12999.7488
16 -19102.4290 -21868.1109
17 -23608.1563 -19102.4290
18 12475.4609 -23608.1563
19 -14366.9931 12475.4609
20 -11020.4817 -14366.9931
21 -26223.6199 -11020.4817
22 -2685.0977 -26223.6199
23 -24948.9035 -2685.0977
24 19136.9475 -24948.9035
25 -34175.6075 19136.9475
26 -17238.4220 -34175.6075
27 5927.3334 -17238.4220
28 999.5841 5927.3334
29 -2218.9730 999.5841
30 10775.3140 -2218.9730
31 6627.6100 10775.3140
32 20184.7555 6627.6100
33 17379.0862 20184.7555
34 8022.1898 17379.0862
35 47249.0959 8022.1898
36 10984.0783 47249.0959
37 -46950.6737 10984.0783
38 11988.9330 -46950.6737
39 12447.0404 11988.9330
40 -4769.7976 12447.0404
41 -12090.3169 -4769.7976
42 10217.9553 -12090.3169
43 19096.4171 10217.9553
44 -24468.3063 19096.4171
45 12856.1676 -24468.3063
46 -40137.8288 12856.1676
47 -20775.7797 -40137.8288
48 7967.6814 -20775.7797
49 81171.5554 7967.6814
50 -21405.0288 81171.5554
51 -11309.2074 -21405.0288
52 3995.8328 -11309.2074
53 25080.4435 3995.8328
54 -50068.4068 25080.4435
55 411.4767 -50068.4068
56 5407.9179 411.4767
57 -1094.3915 5407.9179
58 NA -1094.3915
> dum1 <- dum[2:length(myerror),]
> dum1
lag(myerror, k = 1) myerror
[1,] -21147.7309 -24171.8462
[2,] 13654.7690 -21147.7309
[3,] 14802.9445 13654.7690
[4,] 18876.8096 14802.9445
[5,] 12837.0027 18876.8096
[6,] 16599.6766 12837.0027
[7,] -11768.5107 16599.6766
[8,] 9896.1145 -11768.5107
[9,] -2917.2423 9896.1145
[10,] 34800.7368 -2917.2423
[11,] -1524.4127 34800.7368
[12,] -13916.8610 -1524.4127
[13,] 21102.4566 -13916.8610
[14,] 12999.7488 21102.4566
[15,] -21868.1109 12999.7488
[16,] -19102.4290 -21868.1109
[17,] -23608.1563 -19102.4290
[18,] 12475.4609 -23608.1563
[19,] -14366.9931 12475.4609
[20,] -11020.4817 -14366.9931
[21,] -26223.6199 -11020.4817
[22,] -2685.0977 -26223.6199
[23,] -24948.9035 -2685.0977
[24,] 19136.9475 -24948.9035
[25,] -34175.6075 19136.9475
[26,] -17238.4220 -34175.6075
[27,] 5927.3334 -17238.4220
[28,] 999.5841 5927.3334
[29,] -2218.9730 999.5841
[30,] 10775.3140 -2218.9730
[31,] 6627.6100 10775.3140
[32,] 20184.7555 6627.6100
[33,] 17379.0862 20184.7555
[34,] 8022.1898 17379.0862
[35,] 47249.0959 8022.1898
[36,] 10984.0783 47249.0959
[37,] -46950.6737 10984.0783
[38,] 11988.9330 -46950.6737
[39,] 12447.0404 11988.9330
[40,] -4769.7976 12447.0404
[41,] -12090.3169 -4769.7976
[42,] 10217.9553 -12090.3169
[43,] 19096.4171 10217.9553
[44,] -24468.3063 19096.4171
[45,] 12856.1676 -24468.3063
[46,] -40137.8288 12856.1676
[47,] -20775.7797 -40137.8288
[48,] 7967.6814 -20775.7797
[49,] 81171.5554 7967.6814
[50,] -21405.0288 81171.5554
[51,] -11309.2074 -21405.0288
[52,] 3995.8328 -11309.2074
[53,] 25080.4435 3995.8328
[54,] -50068.4068 25080.4435
[55,] 411.4767 -50068.4068
[56,] 5407.9179 411.4767
[57,] -1094.3915 5407.9179
> z <- as.data.frame(dum1)
> z
lag(myerror, k = 1) myerror
1 -21147.7309 -24171.8462
2 13654.7690 -21147.7309
3 14802.9445 13654.7690
4 18876.8096 14802.9445
5 12837.0027 18876.8096
6 16599.6766 12837.0027
7 -11768.5107 16599.6766
8 9896.1145 -11768.5107
9 -2917.2423 9896.1145
10 34800.7368 -2917.2423
11 -1524.4127 34800.7368
12 -13916.8610 -1524.4127
13 21102.4566 -13916.8610
14 12999.7488 21102.4566
15 -21868.1109 12999.7488
16 -19102.4290 -21868.1109
17 -23608.1563 -19102.4290
18 12475.4609 -23608.1563
19 -14366.9931 12475.4609
20 -11020.4817 -14366.9931
21 -26223.6199 -11020.4817
22 -2685.0977 -26223.6199
23 -24948.9035 -2685.0977
24 19136.9475 -24948.9035
25 -34175.6075 19136.9475
26 -17238.4220 -34175.6075
27 5927.3334 -17238.4220
28 999.5841 5927.3334
29 -2218.9730 999.5841
30 10775.3140 -2218.9730
31 6627.6100 10775.3140
32 20184.7555 6627.6100
33 17379.0862 20184.7555
34 8022.1898 17379.0862
35 47249.0959 8022.1898
36 10984.0783 47249.0959
37 -46950.6737 10984.0783
38 11988.9330 -46950.6737
39 12447.0404 11988.9330
40 -4769.7976 12447.0404
41 -12090.3169 -4769.7976
42 10217.9553 -12090.3169
43 19096.4171 10217.9553
44 -24468.3063 19096.4171
45 12856.1676 -24468.3063
46 -40137.8288 12856.1676
47 -20775.7797 -40137.8288
48 7967.6814 -20775.7797
49 81171.5554 7967.6814
50 -21405.0288 81171.5554
51 -11309.2074 -21405.0288
52 3995.8328 -11309.2074
53 25080.4435 3995.8328
54 -50068.4068 25080.4435
55 411.4767 -50068.4068
56 5407.9179 411.4767
57 -1094.3915 5407.9179
> 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/freestat/rcomp/tmp/7gr9z1292511296.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/freestat/rcomp/tmp/89irk1292511296.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/freestat/rcomp/tmp/99irk1292511296.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/freestat/rcomp/tmp/109irk1292511296.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/freestat/rcomp/createtable file can be downloaded at http://www.wessa.net/cretab
> load(file="/var/www/html/freestat/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/freestat/rcomp/tmp/11c07p1292511296.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/freestat/rcomp/tmp/12y16d1292511296.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/freestat/rcomp/tmp/1352l71292511296.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/freestat/rcomp/tmp/14xbka1292511296.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/freestat/rcomp/tmp/150u0g1292511296.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/freestat/rcomp/tmp/16fly71292511296.tab")
+ }
>
> try(system("convert tmp/1v8bt1292511296.ps tmp/1v8bt1292511296.png",intern=TRUE))
character(0)
> try(system("convert tmp/2v8bt1292511296.ps tmp/2v8bt1292511296.png",intern=TRUE))
character(0)
> try(system("convert tmp/3v8bt1292511296.ps tmp/3v8bt1292511296.png",intern=TRUE))
character(0)
> try(system("convert tmp/46hse1292511296.ps tmp/46hse1292511296.png",intern=TRUE))
character(0)
> try(system("convert tmp/56hse1292511296.ps tmp/56hse1292511296.png",intern=TRUE))
character(0)
> try(system("convert tmp/66hse1292511296.ps tmp/66hse1292511296.png",intern=TRUE))
character(0)
> try(system("convert tmp/7gr9z1292511296.ps tmp/7gr9z1292511296.png",intern=TRUE))
character(0)
> try(system("convert tmp/89irk1292511296.ps tmp/89irk1292511296.png",intern=TRUE))
character(0)
> try(system("convert tmp/99irk1292511296.ps tmp/99irk1292511296.png",intern=TRUE))
character(0)
> try(system("convert tmp/109irk1292511296.ps tmp/109irk1292511296.png",intern=TRUE))
character(0)
>
>
> proc.time()
user system elapsed
3.855 2.544 4.734