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(989236.00
+ ,10489.94
+ ,1008380.00
+ ,10766.23
+ ,1207763.00
+ ,10503.76
+ ,1368839.00
+ ,10192.51
+ ,1469798.00
+ ,10467.48
+ ,1498721.00
+ ,10274.97
+ ,1761769.00
+ ,10640.91
+ ,1653214.00
+ ,10481.60
+ ,1599104.00
+ ,10568.70
+ ,1421179.00
+ ,10440.07
+ ,1163995.00
+ ,10805.87
+ ,1037735.00
+ ,10717.50
+ ,1015407.00
+ ,10864.86
+ ,1039210.00
+ ,10993.41
+ ,1258049.00
+ ,11109.32
+ ,1469445.00
+ ,11367.14
+ ,1552346.00
+ ,11168.31
+ ,1549144.00
+ ,11150.22
+ ,1785895.00
+ ,11185.68
+ ,1662335.00
+ ,11381.15
+ ,1629440.00
+ ,11679.07
+ ,1467430.00
+ ,12080.73
+ ,1202209.00
+ ,12221.93
+ ,1076982.00
+ ,12463.15
+ ,1039367.00
+ ,12621.69
+ ,1063449.00
+ ,12268.63
+ ,1335135.00
+ ,12354.35
+ ,1491602.00
+ ,13062.91
+ ,1591972.00
+ ,13627.64
+ ,1641248.00
+ ,13408.62
+ ,1898849.00
+ ,13211.99
+ ,1798580.00
+ ,13357.74
+ ,1762444.00
+ ,13895.63
+ ,1622044.00
+ ,13930.01
+ ,1368955.00
+ ,13371.72
+ ,1262973.00
+ ,13264.82
+ ,1195650.00
+ ,12650.36
+ ,1269530.00
+ ,12266.39
+ ,1479279.00
+ ,12262.89
+ ,1607819.00
+ ,12820.13
+ ,1712466.00
+ ,12638.32
+ ,1721766.00
+ ,11350.01
+ ,1949843.00
+ ,11378.02
+ ,1821326.00
+ ,11543.55
+ ,1757802.00
+ ,10850.66
+ ,1590367.00
+ ,9325.01
+ ,1260647.00
+ ,8829.04
+ ,1149235.00
+ ,8776.39
+ ,1016367.00
+ ,8000.86
+ ,1027885.00
+ ,7062.93
+ ,1262159.00
+ ,7608.92
+ ,1520854.00
+ ,8168.12
+ ,1544144.00
+ ,8500.33
+ ,1564709.00
+ ,8447.00
+ ,1821776.00
+ ,9171.61
+ ,1741365.00
+ ,9496.28
+ ,1623386.00
+ ,9712.28
+ ,1498658.00
+ ,9712.73
+ ,1241822.00
+ ,10344.84
+ ,1136029.00
+ ,10428.05)
+ ,dim=c(2
+ ,60)
+ ,dimnames=list(c('Passengersbrussels'
+ ,'DJIA')
+ ,1:60))
> y <- array(NA,dim=c(2,60),dimnames=list(c('Passengersbrussels','DJIA'),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 = '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
Passengersbrussels DJIA M1 M2 M3 M4 M5 M6 M7 M8 M9 M10 M11 t
1 989236 10489.94 1 0 0 0 0 0 0 0 0 0 0 1
2 1008380 10766.23 0 1 0 0 0 0 0 0 0 0 0 2
3 1207763 10503.76 0 0 1 0 0 0 0 0 0 0 0 3
4 1368839 10192.51 0 0 0 1 0 0 0 0 0 0 0 4
5 1469798 10467.48 0 0 0 0 1 0 0 0 0 0 0 5
6 1498721 10274.97 0 0 0 0 0 1 0 0 0 0 0 6
7 1761769 10640.91 0 0 0 0 0 0 1 0 0 0 0 7
8 1653214 10481.60 0 0 0 0 0 0 0 1 0 0 0 8
9 1599104 10568.70 0 0 0 0 0 0 0 0 1 0 0 9
10 1421179 10440.07 0 0 0 0 0 0 0 0 0 1 0 10
11 1163995 10805.87 0 0 0 0 0 0 0 0 0 0 1 11
12 1037735 10717.50 0 0 0 0 0 0 0 0 0 0 0 12
13 1015407 10864.86 1 0 0 0 0 0 0 0 0 0 0 13
14 1039210 10993.41 0 1 0 0 0 0 0 0 0 0 0 14
15 1258049 11109.32 0 0 1 0 0 0 0 0 0 0 0 15
16 1469445 11367.14 0 0 0 1 0 0 0 0 0 0 0 16
17 1552346 11168.31 0 0 0 0 1 0 0 0 0 0 0 17
18 1549144 11150.22 0 0 0 0 0 1 0 0 0 0 0 18
19 1785895 11185.68 0 0 0 0 0 0 1 0 0 0 0 19
20 1662335 11381.15 0 0 0 0 0 0 0 1 0 0 0 20
21 1629440 11679.07 0 0 0 0 0 0 0 0 1 0 0 21
22 1467430 12080.73 0 0 0 0 0 0 0 0 0 1 0 22
23 1202209 12221.93 0 0 0 0 0 0 0 0 0 0 1 23
24 1076982 12463.15 0 0 0 0 0 0 0 0 0 0 0 24
25 1039367 12621.69 1 0 0 0 0 0 0 0 0 0 0 25
26 1063449 12268.63 0 1 0 0 0 0 0 0 0 0 0 26
27 1335135 12354.35 0 0 1 0 0 0 0 0 0 0 0 27
28 1491602 13062.91 0 0 0 1 0 0 0 0 0 0 0 28
29 1591972 13627.64 0 0 0 0 1 0 0 0 0 0 0 29
30 1641248 13408.62 0 0 0 0 0 1 0 0 0 0 0 30
31 1898849 13211.99 0 0 0 0 0 0 1 0 0 0 0 31
32 1798580 13357.74 0 0 0 0 0 0 0 1 0 0 0 32
33 1762444 13895.63 0 0 0 0 0 0 0 0 1 0 0 33
34 1622044 13930.01 0 0 0 0 0 0 0 0 0 1 0 34
35 1368955 13371.72 0 0 0 0 0 0 0 0 0 0 1 35
36 1262973 13264.82 0 0 0 0 0 0 0 0 0 0 0 36
37 1195650 12650.36 1 0 0 0 0 0 0 0 0 0 0 37
38 1269530 12266.39 0 1 0 0 0 0 0 0 0 0 0 38
39 1479279 12262.89 0 0 1 0 0 0 0 0 0 0 0 39
40 1607819 12820.13 0 0 0 1 0 0 0 0 0 0 0 40
41 1712466 12638.32 0 0 0 0 1 0 0 0 0 0 0 41
42 1721766 11350.01 0 0 0 0 0 1 0 0 0 0 0 42
43 1949843 11378.02 0 0 0 0 0 0 1 0 0 0 0 43
44 1821326 11543.55 0 0 0 0 0 0 0 1 0 0 0 44
45 1757802 10850.66 0 0 0 0 0 0 0 0 1 0 0 45
46 1590367 9325.01 0 0 0 0 0 0 0 0 0 1 0 46
47 1260647 8829.04 0 0 0 0 0 0 0 0 0 0 1 47
48 1149235 8776.39 0 0 0 0 0 0 0 0 0 0 0 48
49 1016367 8000.86 1 0 0 0 0 0 0 0 0 0 0 49
50 1027885 7062.93 0 1 0 0 0 0 0 0 0 0 0 50
51 1262159 7608.92 0 0 1 0 0 0 0 0 0 0 0 51
52 1520854 8168.12 0 0 0 1 0 0 0 0 0 0 0 52
53 1544144 8500.33 0 0 0 0 1 0 0 0 0 0 0 53
54 1564709 8447.00 0 0 0 0 0 1 0 0 0 0 0 54
55 1821776 9171.61 0 0 0 0 0 0 1 0 0 0 0 55
56 1741365 9496.28 0 0 0 0 0 0 0 1 0 0 0 56
57 1623386 9712.28 0 0 0 0 0 0 0 0 1 0 0 57
58 1498658 9712.73 0 0 0 0 0 0 0 0 0 1 0 58
59 1241822 10344.84 0 0 0 0 0 0 0 0 0 0 1 59
60 1136029 10428.05 0 0 0 0 0 0 0 0 0 0 0 60
> k <- length(x[1,])
> df <- as.data.frame(x)
> (mylm <- lm(df))
Call:
lm(formula = df)
Coefficients:
(Intercept) DJIA M1 M2 M3 M4
643506.89 33.06 -37601.18 -2084.93 218151.00 386307.76
M5 M6 M7 M8 M9 M10
460143.96 489460.30 728273.67 612201.85 544958.34 395143.63
M11 t
118806.65 3366.03
> (mysum <- summary(mylm))
Call:
lm(formula = df)
Residuals:
Min 1Q Median 3Q Max
-77987 -26649 -3412 31547 94728
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 643506.890 55597.882 11.574 3.19e-15 ***
DJIA 33.055 4.053 8.155 1.75e-10 ***
M1 -37601.179 31091.216 -1.209 0.232695
M2 -2084.931 31121.009 -0.067 0.946877
M3 218151.003 31033.333 7.030 8.23e-09 ***
M4 386307.757 30907.224 12.499 < 2e-16 ***
M5 460143.964 30859.092 14.911 < 2e-16 ***
M6 489460.296 30867.966 15.857 < 2e-16 ***
M7 728273.669 30811.808 23.636 < 2e-16 ***
M8 612201.854 30785.045 19.886 < 2e-16 ***
M9 544958.342 30773.524 17.709 < 2e-16 ***
M10 395143.632 30760.584 12.846 < 2e-16 ***
M11 118806.648 30752.229 3.863 0.000348 ***
t 3366.033 388.608 8.662 3.19e-11 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 48620 on 46 degrees of freedom
Multiple R-squared: 0.9746, Adjusted R-squared: 0.9675
F-statistic: 136 on 13 and 46 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,] 5.877554e-03 0.0117551072 0.99412245
[2,] 3.570381e-03 0.0071407613 0.99642962
[3,] 1.416502e-03 0.0028330050 0.99858350
[4,] 5.076172e-03 0.0101523447 0.99492383
[5,] 2.842925e-03 0.0056858500 0.99715708
[6,] 1.117094e-03 0.0022341884 0.99888291
[7,] 3.835236e-04 0.0007670471 0.99961648
[8,] 1.480505e-04 0.0002961010 0.99985195
[9,] 1.098085e-04 0.0002196171 0.99989019
[10,] 8.623863e-05 0.0001724773 0.99991376
[11,] 1.397717e-04 0.0002795434 0.99986023
[12,] 1.427064e-04 0.0002854129 0.99985729
[13,] 1.698663e-04 0.0003397326 0.99983013
[14,] 5.543466e-04 0.0011086933 0.99944565
[15,] 1.807458e-03 0.0036149161 0.99819254
[16,] 9.235501e-03 0.0184710022 0.99076450
[17,] 1.894380e-02 0.0378876054 0.98105620
[18,] 1.366452e-01 0.2732903255 0.86335484
[19,] 3.908119e-01 0.7816237703 0.60918811
[20,] 7.368739e-01 0.5262521348 0.26312607
[21,] 7.022275e-01 0.5955450771 0.29777254
[22,] 7.922113e-01 0.4155774854 0.20778874
[23,] 8.184650e-01 0.3630699927 0.18153500
[24,] 9.539826e-01 0.0920347059 0.04601735
[25,] 9.036036e-01 0.1927928245 0.09639641
[26,] 8.170303e-01 0.3659394035 0.18296970
[27,] 6.718327e-01 0.6563346125 0.32816731
> postscript(file="/var/www/html/rcomp/tmp/18z9j1292107217.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/2iqq41292107217.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/3iqq41292107217.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/4iqq41292107217.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/5biqp1292107217.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 = 60
Frequency = 1
1 2 3 4 5 6
33214.4153 4343.2396 -11199.6597 -11357.9324 3309.5660 5913.7090
7 8 9 10 11 12
14685.9865 24102.8343 30991.1830 3766.7846 7462.0468 -436.2271
13 14 15 16 17 18
6599.8661 -12728.6953 -21323.1213 10027.7327 22298.9103 -12987.4808
19 20 21 22 23 24
-19588.0333 -36903.6034 -15769.0077 -44607.3886 -41524.8689 -59284.8934
25 26 27 28 29 30
-67905.3604 -71035.0794 -25784.5609 -64262.1271 -59761.7801 -35928.3367
31 32 33 34 35 36
-14007.0464 -6388.0987 3573.1770 8485.4073 46821.8939 59814.1378
37 38 39 40 41 42
47037.5451 94727.5705 80990.2975 19587.6843 53042.2579 72245.5791
43 44 45 46 47 48
57217.2898 35934.4004 59191.6811 88636.4279 48281.8980 54050.8830
49 50 51 52 53 54
-18946.4661 -15307.0354 -22682.9557 46004.6424 -18888.9540 -29243.4706
55 56 57 58 59 60
-38308.1966 -16745.5327 -77987.0334 -56281.2312 -61040.9698 -54143.9003
> postscript(file="/var/www/html/rcomp/tmp/6biqp1292107217.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 = 60
Frequency = 1
lag(myerror, k = 1) myerror
0 33214.4153 NA
1 4343.2396 33214.4153
2 -11199.6597 4343.2396
3 -11357.9324 -11199.6597
4 3309.5660 -11357.9324
5 5913.7090 3309.5660
6 14685.9865 5913.7090
7 24102.8343 14685.9865
8 30991.1830 24102.8343
9 3766.7846 30991.1830
10 7462.0468 3766.7846
11 -436.2271 7462.0468
12 6599.8661 -436.2271
13 -12728.6953 6599.8661
14 -21323.1213 -12728.6953
15 10027.7327 -21323.1213
16 22298.9103 10027.7327
17 -12987.4808 22298.9103
18 -19588.0333 -12987.4808
19 -36903.6034 -19588.0333
20 -15769.0077 -36903.6034
21 -44607.3886 -15769.0077
22 -41524.8689 -44607.3886
23 -59284.8934 -41524.8689
24 -67905.3604 -59284.8934
25 -71035.0794 -67905.3604
26 -25784.5609 -71035.0794
27 -64262.1271 -25784.5609
28 -59761.7801 -64262.1271
29 -35928.3367 -59761.7801
30 -14007.0464 -35928.3367
31 -6388.0987 -14007.0464
32 3573.1770 -6388.0987
33 8485.4073 3573.1770
34 46821.8939 8485.4073
35 59814.1378 46821.8939
36 47037.5451 59814.1378
37 94727.5705 47037.5451
38 80990.2975 94727.5705
39 19587.6843 80990.2975
40 53042.2579 19587.6843
41 72245.5791 53042.2579
42 57217.2898 72245.5791
43 35934.4004 57217.2898
44 59191.6811 35934.4004
45 88636.4279 59191.6811
46 48281.8980 88636.4279
47 54050.8830 48281.8980
48 -18946.4661 54050.8830
49 -15307.0354 -18946.4661
50 -22682.9557 -15307.0354
51 46004.6424 -22682.9557
52 -18888.9540 46004.6424
53 -29243.4706 -18888.9540
54 -38308.1966 -29243.4706
55 -16745.5327 -38308.1966
56 -77987.0334 -16745.5327
57 -56281.2312 -77987.0334
58 -61040.9698 -56281.2312
59 -54143.9003 -61040.9698
60 NA -54143.9003
> dum1 <- dum[2:length(myerror),]
> dum1
lag(myerror, k = 1) myerror
[1,] 4343.2396 33214.4153
[2,] -11199.6597 4343.2396
[3,] -11357.9324 -11199.6597
[4,] 3309.5660 -11357.9324
[5,] 5913.7090 3309.5660
[6,] 14685.9865 5913.7090
[7,] 24102.8343 14685.9865
[8,] 30991.1830 24102.8343
[9,] 3766.7846 30991.1830
[10,] 7462.0468 3766.7846
[11,] -436.2271 7462.0468
[12,] 6599.8661 -436.2271
[13,] -12728.6953 6599.8661
[14,] -21323.1213 -12728.6953
[15,] 10027.7327 -21323.1213
[16,] 22298.9103 10027.7327
[17,] -12987.4808 22298.9103
[18,] -19588.0333 -12987.4808
[19,] -36903.6034 -19588.0333
[20,] -15769.0077 -36903.6034
[21,] -44607.3886 -15769.0077
[22,] -41524.8689 -44607.3886
[23,] -59284.8934 -41524.8689
[24,] -67905.3604 -59284.8934
[25,] -71035.0794 -67905.3604
[26,] -25784.5609 -71035.0794
[27,] -64262.1271 -25784.5609
[28,] -59761.7801 -64262.1271
[29,] -35928.3367 -59761.7801
[30,] -14007.0464 -35928.3367
[31,] -6388.0987 -14007.0464
[32,] 3573.1770 -6388.0987
[33,] 8485.4073 3573.1770
[34,] 46821.8939 8485.4073
[35,] 59814.1378 46821.8939
[36,] 47037.5451 59814.1378
[37,] 94727.5705 47037.5451
[38,] 80990.2975 94727.5705
[39,] 19587.6843 80990.2975
[40,] 53042.2579 19587.6843
[41,] 72245.5791 53042.2579
[42,] 57217.2898 72245.5791
[43,] 35934.4004 57217.2898
[44,] 59191.6811 35934.4004
[45,] 88636.4279 59191.6811
[46,] 48281.8980 88636.4279
[47,] 54050.8830 48281.8980
[48,] -18946.4661 54050.8830
[49,] -15307.0354 -18946.4661
[50,] -22682.9557 -15307.0354
[51,] 46004.6424 -22682.9557
[52,] -18888.9540 46004.6424
[53,] -29243.4706 -18888.9540
[54,] -38308.1966 -29243.4706
[55,] -16745.5327 -38308.1966
[56,] -77987.0334 -16745.5327
[57,] -56281.2312 -77987.0334
[58,] -61040.9698 -56281.2312
[59,] -54143.9003 -61040.9698
> z <- as.data.frame(dum1)
> z
lag(myerror, k = 1) myerror
1 4343.2396 33214.4153
2 -11199.6597 4343.2396
3 -11357.9324 -11199.6597
4 3309.5660 -11357.9324
5 5913.7090 3309.5660
6 14685.9865 5913.7090
7 24102.8343 14685.9865
8 30991.1830 24102.8343
9 3766.7846 30991.1830
10 7462.0468 3766.7846
11 -436.2271 7462.0468
12 6599.8661 -436.2271
13 -12728.6953 6599.8661
14 -21323.1213 -12728.6953
15 10027.7327 -21323.1213
16 22298.9103 10027.7327
17 -12987.4808 22298.9103
18 -19588.0333 -12987.4808
19 -36903.6034 -19588.0333
20 -15769.0077 -36903.6034
21 -44607.3886 -15769.0077
22 -41524.8689 -44607.3886
23 -59284.8934 -41524.8689
24 -67905.3604 -59284.8934
25 -71035.0794 -67905.3604
26 -25784.5609 -71035.0794
27 -64262.1271 -25784.5609
28 -59761.7801 -64262.1271
29 -35928.3367 -59761.7801
30 -14007.0464 -35928.3367
31 -6388.0987 -14007.0464
32 3573.1770 -6388.0987
33 8485.4073 3573.1770
34 46821.8939 8485.4073
35 59814.1378 46821.8939
36 47037.5451 59814.1378
37 94727.5705 47037.5451
38 80990.2975 94727.5705
39 19587.6843 80990.2975
40 53042.2579 19587.6843
41 72245.5791 53042.2579
42 57217.2898 72245.5791
43 35934.4004 57217.2898
44 59191.6811 35934.4004
45 88636.4279 59191.6811
46 48281.8980 88636.4279
47 54050.8830 48281.8980
48 -18946.4661 54050.8830
49 -15307.0354 -18946.4661
50 -22682.9557 -15307.0354
51 46004.6424 -22682.9557
52 -18888.9540 46004.6424
53 -29243.4706 -18888.9540
54 -38308.1966 -29243.4706
55 -16745.5327 -38308.1966
56 -77987.0334 -16745.5327
57 -56281.2312 -77987.0334
58 -61040.9698 -56281.2312
59 -54143.9003 -61040.9698
> 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/74rpa1292107217.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/84rpa1292107217.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/9x06v1292107217.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/10x06v1292107217.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/11iin11292107217.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/12l13p1292107217.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/13s2ii1292107217.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/14lbhl1292107217.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/15ocy91292107217.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/162mwi1292107217.tab")
+ }
>
> try(system("convert tmp/18z9j1292107217.ps tmp/18z9j1292107217.png",intern=TRUE))
character(0)
> try(system("convert tmp/2iqq41292107217.ps tmp/2iqq41292107217.png",intern=TRUE))
character(0)
> try(system("convert tmp/3iqq41292107217.ps tmp/3iqq41292107217.png",intern=TRUE))
character(0)
> try(system("convert tmp/4iqq41292107217.ps tmp/4iqq41292107217.png",intern=TRUE))
character(0)
> try(system("convert tmp/5biqp1292107217.ps tmp/5biqp1292107217.png",intern=TRUE))
character(0)
> try(system("convert tmp/6biqp1292107217.ps tmp/6biqp1292107217.png",intern=TRUE))
character(0)
> try(system("convert tmp/74rpa1292107217.ps tmp/74rpa1292107217.png",intern=TRUE))
character(0)
> try(system("convert tmp/84rpa1292107217.ps tmp/84rpa1292107217.png",intern=TRUE))
character(0)
> try(system("convert tmp/9x06v1292107217.ps tmp/9x06v1292107217.png",intern=TRUE))
character(0)
> try(system("convert tmp/10x06v1292107217.ps tmp/10x06v1292107217.png",intern=TRUE))
character(0)
>
>
> proc.time()
user system elapsed
2.437 1.634 11.093