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(267413,21.4,267366,26.4,264777,26.4,258863,29.4,254844,34.4,254868,24.4,277267,26.4,285351,25.4,286602,31.4,283042,27.4,276687,27.4,277915,29.4,277128,32.4,277103,26.4,275037,22.4,270150,19.4,267140,21.4,264993,23.4,287259,23.4,291186,25.4,292300,28.4,288186,27.4,281477,21.4,282656,17.4,280190,24.4,280408,26.4,276836,22.4,275216,14.4,274352,18.4,271311,25.4,289802,29.4,290726,26.4,292300,26.4,278506,20.4,269826,26.4,265861,29.4,269034,33.4,264176,32.4,255198,35.4,253353,34.4,246057,36.4,235372,32.4,258556,34.4,260993,31.4,254663,27.4,250643,27.4,243422,30.4,247105,32.4,248541,32.4,245039,27.4,237080,31.4,237085,29.4,225554,27.4,226839,25.4,247934,26.4,248333,23.4,246969,18.4,245098,22.4,246263,17.4,255765,17.4,264319,11.4,268347,9.4,273046,6.4,273963,0,267430,7.8,271993,7.9,292710,12,295881,16.9,293299,12.3),dim=c(2,69),dimnames=list(c('Y','X'),1:69))
> y <- array(NA,dim=c(2,69),dimnames=list(c('Y','X'),1:69))
> 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 X M1 M2 M3 M4 M5 M6 M7 M8 M9 M10 M11 t
1 267413 21.4 1 0 0 0 0 0 0 0 0 0 0 1
2 267366 26.4 0 1 0 0 0 0 0 0 0 0 0 2
3 264777 26.4 0 0 1 0 0 0 0 0 0 0 0 3
4 258863 29.4 0 0 0 1 0 0 0 0 0 0 0 4
5 254844 34.4 0 0 0 0 1 0 0 0 0 0 0 5
6 254868 24.4 0 0 0 0 0 1 0 0 0 0 0 6
7 277267 26.4 0 0 0 0 0 0 1 0 0 0 0 7
8 285351 25.4 0 0 0 0 0 0 0 1 0 0 0 8
9 286602 31.4 0 0 0 0 0 0 0 0 1 0 0 9
10 283042 27.4 0 0 0 0 0 0 0 0 0 1 0 10
11 276687 27.4 0 0 0 0 0 0 0 0 0 0 1 11
12 277915 29.4 0 0 0 0 0 0 0 0 0 0 0 12
13 277128 32.4 1 0 0 0 0 0 0 0 0 0 0 13
14 277103 26.4 0 1 0 0 0 0 0 0 0 0 0 14
15 275037 22.4 0 0 1 0 0 0 0 0 0 0 0 15
16 270150 19.4 0 0 0 1 0 0 0 0 0 0 0 16
17 267140 21.4 0 0 0 0 1 0 0 0 0 0 0 17
18 264993 23.4 0 0 0 0 0 1 0 0 0 0 0 18
19 287259 23.4 0 0 0 0 0 0 1 0 0 0 0 19
20 291186 25.4 0 0 0 0 0 0 0 1 0 0 0 20
21 292300 28.4 0 0 0 0 0 0 0 0 1 0 0 21
22 288186 27.4 0 0 0 0 0 0 0 0 0 1 0 22
23 281477 21.4 0 0 0 0 0 0 0 0 0 0 1 23
24 282656 17.4 0 0 0 0 0 0 0 0 0 0 0 24
25 280190 24.4 1 0 0 0 0 0 0 0 0 0 0 25
26 280408 26.4 0 1 0 0 0 0 0 0 0 0 0 26
27 276836 22.4 0 0 1 0 0 0 0 0 0 0 0 27
28 275216 14.4 0 0 0 1 0 0 0 0 0 0 0 28
29 274352 18.4 0 0 0 0 1 0 0 0 0 0 0 29
30 271311 25.4 0 0 0 0 0 1 0 0 0 0 0 30
31 289802 29.4 0 0 0 0 0 0 1 0 0 0 0 31
32 290726 26.4 0 0 0 0 0 0 0 1 0 0 0 32
33 292300 26.4 0 0 0 0 0 0 0 0 1 0 0 33
34 278506 20.4 0 0 0 0 0 0 0 0 0 1 0 34
35 269826 26.4 0 0 0 0 0 0 0 0 0 0 1 35
36 265861 29.4 0 0 0 0 0 0 0 0 0 0 0 36
37 269034 33.4 1 0 0 0 0 0 0 0 0 0 0 37
38 264176 32.4 0 1 0 0 0 0 0 0 0 0 0 38
39 255198 35.4 0 0 1 0 0 0 0 0 0 0 0 39
40 253353 34.4 0 0 0 1 0 0 0 0 0 0 0 40
41 246057 36.4 0 0 0 0 1 0 0 0 0 0 0 41
42 235372 32.4 0 0 0 0 0 1 0 0 0 0 0 42
43 258556 34.4 0 0 0 0 0 0 1 0 0 0 0 43
44 260993 31.4 0 0 0 0 0 0 0 1 0 0 0 44
45 254663 27.4 0 0 0 0 0 0 0 0 1 0 0 45
46 250643 27.4 0 0 0 0 0 0 0 0 0 1 0 46
47 243422 30.4 0 0 0 0 0 0 0 0 0 0 1 47
48 247105 32.4 0 0 0 0 0 0 0 0 0 0 0 48
49 248541 32.4 1 0 0 0 0 0 0 0 0 0 0 49
50 245039 27.4 0 1 0 0 0 0 0 0 0 0 0 50
51 237080 31.4 0 0 1 0 0 0 0 0 0 0 0 51
52 237085 29.4 0 0 0 1 0 0 0 0 0 0 0 52
53 225554 27.4 0 0 0 0 1 0 0 0 0 0 0 53
54 226839 25.4 0 0 0 0 0 1 0 0 0 0 0 54
55 247934 26.4 0 0 0 0 0 0 1 0 0 0 0 55
56 248333 23.4 0 0 0 0 0 0 0 1 0 0 0 56
57 246969 18.4 0 0 0 0 0 0 0 0 1 0 0 57
58 245098 22.4 0 0 0 0 0 0 0 0 0 1 0 58
59 246263 17.4 0 0 0 0 0 0 0 0 0 0 1 59
60 255765 17.4 0 0 0 0 0 0 0 0 0 0 0 60
61 264319 11.4 1 0 0 0 0 0 0 0 0 0 0 61
62 268347 9.4 0 1 0 0 0 0 0 0 0 0 0 62
63 273046 6.4 0 0 1 0 0 0 0 0 0 0 0 63
64 273963 0.0 0 0 0 1 0 0 0 0 0 0 0 64
65 267430 7.8 0 0 0 0 1 0 0 0 0 0 0 65
66 271993 7.9 0 0 0 0 0 1 0 0 0 0 0 66
67 292710 12.0 0 0 0 0 0 0 1 0 0 0 0 67
68 295881 16.9 0 0 0 0 0 0 0 1 0 0 0 68
69 293299 12.3 0 0 0 0 0 0 0 0 1 0 0 69
> k <- length(x[1,])
> df <- as.data.frame(x)
> (mylm <- lm(df))
Call:
lm(formula = df)
Coefficients:
(Intercept) X M1 M2 M3 M4
318115.4 -1319.5 194.8 -1514.4 -5277.1 -10799.9
M5 M6 M7 M8 M9 M10
-11679.7 -14336.1 10431.4 13434.5 11894.5 1915.0
M11 t
-3645.0 -527.9
> (mysum <- summary(mylm))
Call:
lm(formula = df)
Residuals:
Min 1Q Median 3Q Max
-28674 -7547 430 8433 22525
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 318115.40 9097.36 34.968 < 2e-16 ***
X -1319.54 216.99 -6.081 1.19e-07 ***
M1 194.85 7453.37 0.026 0.9792
M2 -1514.42 7453.87 -0.203 0.8397
M3 -5277.10 7455.59 -0.708 0.4821
M4 -10799.90 7505.21 -1.439 0.1558
M5 -11679.67 7447.88 -1.568 0.1226
M6 -14336.12 7457.22 -1.922 0.0597 .
M7 10431.39 7444.60 1.401 0.1668
M8 13434.48 7445.46 1.804 0.0766 .
M9 11894.52 7448.73 1.597 0.1160
M10 1914.99 7777.21 0.246 0.8064
M11 -3644.97 7777.06 -0.469 0.6411
t -527.85 82.61 -6.390 3.75e-08 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 12290 on 55 degrees of freedom
Multiple R-squared: 0.5963, Adjusted R-squared: 0.5009
F-statistic: 6.25 on 13 and 55 DF, p-value: 5.139e-07
> 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,] 7.478439e-05 1.495688e-04 0.9999252156
[2,] 2.634064e-06 5.268128e-06 0.9999973659
[3,] 1.350711e-07 2.701422e-07 0.9999998649
[4,] 5.209540e-06 1.041908e-05 0.9999947905
[5,] 4.282942e-06 8.565885e-06 0.9999957171
[6,] 2.065935e-06 4.131870e-06 0.9999979341
[7,] 1.209564e-06 2.419127e-06 0.9999987904
[8,] 5.986398e-07 1.197280e-06 0.9999994014
[9,] 2.963527e-07 5.927053e-07 0.9999997036
[10,] 8.104860e-08 1.620972e-07 0.9999999190
[11,] 3.097778e-08 6.195556e-08 0.9999999690
[12,] 7.124200e-09 1.424840e-08 0.9999999929
[13,] 1.845534e-09 3.691069e-09 0.9999999982
[14,] 3.209554e-10 6.419109e-10 0.9999999997
[15,] 6.800887e-11 1.360177e-10 0.9999999999
[16,] 2.224326e-10 4.448652e-10 0.9999999998
[17,] 4.632970e-10 9.265940e-10 0.9999999995
[18,] 1.538922e-07 3.077845e-07 0.9999998461
[19,] 2.386658e-06 4.773317e-06 0.9999976133
[20,] 1.820185e-05 3.640369e-05 0.9999817982
[21,] 1.905714e-05 3.811428e-05 0.9999809429
[22,] 3.814777e-05 7.629554e-05 0.9999618522
[23,] 6.237067e-05 1.247413e-04 0.9999376293
[24,] 7.338233e-05 1.467647e-04 0.9999266177
[25,] 2.182518e-04 4.365037e-04 0.9997817482
[26,] 1.040192e-03 2.080384e-03 0.9989598081
[27,] 2.186648e-03 4.373295e-03 0.9978133523
[28,] 5.537757e-03 1.107551e-02 0.9944622428
[29,] 4.433513e-02 8.867025e-02 0.9556648749
[30,] 1.959352e-01 3.918704e-01 0.8040648010
[31,] 4.854598e-01 9.709196e-01 0.5145401983
[32,] 8.133929e-01 3.732142e-01 0.1866070925
[33,] 9.501974e-01 9.960529e-02 0.0498026473
[34,] 9.998668e-01 2.663993e-04 0.0001331996
[35,] 9.998093e-01 3.813191e-04 0.0001906596
[36,] 9.991861e-01 1.627722e-03 0.0008138611
> postscript(file="/var/www/html/rcomp/tmp/13vcg1258715003.ps",horizontal=F,pagecentre=F,paper="special",width=8.3333333333333,height=5.5555555555556)
> plot(x[,1], type='l', main='Actuals and Interpolation', ylab='value of Actuals and Interpolation (dots)', xlab='time or index')
> points(x[,1]-mysum$resid)
> grid()
> dev.off()
null device
1
> postscript(file="/var/www/html/rcomp/tmp/2iehl1258715003.ps",horizontal=F,pagecentre=F,paper="special",width=8.3333333333333,height=5.5555555555556)
> plot(mysum$resid, type='b', pch=19, main='Residuals', ylab='value of Residuals', xlab='time or index')
> grid()
> dev.off()
null device
1
> postscript(file="/var/www/html/rcomp/tmp/3leyy1258715003.ps",horizontal=F,pagecentre=F,paper="special",width=8.3333333333333,height=5.5555555555556)
> hist(mysum$resid, main='Residual Histogram', xlab='values of Residuals')
> grid()
> dev.off()
null device
1
> postscript(file="/var/www/html/rcomp/tmp/4046s1258715003.ps",horizontal=F,pagecentre=F,paper="special",width=8.3333333333333,height=5.5555555555556)
> densityplot(~mysum$resid,col='black',main='Residual Density Plot', xlab='values of Residuals')
> dev.off()
null device
1
> postscript(file="/var/www/html/rcomp/tmp/5f02p1258715003.ps",horizontal=F,pagecentre=F,paper="special",width=8.3333333333333,height=5.5555555555556)
> qqnorm(mysum$resid, main='Residual Normal Q-Q Plot')
> qqline(mysum$resid)
> grid()
> dev.off()
null device
1
> (myerror <- as.ts(mysum$resid))
Time Series:
Start = 1
End = 69
Frequency = 1
1 2 3 4 5 6
-22131.3064 -13343.4955 -11641.9709 -7546.7033 -3560.4011 -13547.4692
7 8 9 10 11 12
-12749.0510 -8459.8272 2776.2058 4445.4429 4178.2577 4928.2093
13 14 15 16 17 18
8432.8223 2727.7273 -325.8957 -3120.8494 -2084.1579 1592.2167
19 20 21 22 23 24
-381.4389 3709.3956 10849.8179 15923.6657 7385.2591 168.9894
25 26 27 28 29 30
7272.7499 12366.9501 7807.3271 1681.6889 7503.4543 16883.5133
31 32 33 34 35 36
16413.0053 10903.1553 14544.9670 3341.1303 8666.1664 5542.6550
37 38 39 40 41 42
14326.8048 10386.3943 9657.5295 12543.6496 9294.3412 -3484.5057
43 44 45 46 47 48
-1901.0874 -5897.9374 -15438.2733 -8950.8886 -6125.4632 -2920.5116
49 50 51 52 53 54
-1151.5093 -9014.0674 -7404.3952 -3987.8120 -16750.2681 -14920.0411
55 56 57 58 59 60
-16745.1598 -22780.0097 -28673.8826 -14759.3503 -14104.2200 -7719.3422
61 62 63 64 65 66
-6749.5613 -3123.5087 1907.4052 430.0261 5597.0316 13476.2861
67 68 69
15363.7318 22525.2233 15941.1652
> postscript(file="/var/www/html/rcomp/tmp/6j5b91258715003.ps",horizontal=F,pagecentre=F,paper="special",width=8.3333333333333,height=5.5555555555556)
> dum <- cbind(lag(myerror,k=1),myerror)
> dum
Time Series:
Start = 0
End = 69
Frequency = 1
lag(myerror, k = 1) myerror
0 -22131.3064 NA
1 -13343.4955 -22131.3064
2 -11641.9709 -13343.4955
3 -7546.7033 -11641.9709
4 -3560.4011 -7546.7033
5 -13547.4692 -3560.4011
6 -12749.0510 -13547.4692
7 -8459.8272 -12749.0510
8 2776.2058 -8459.8272
9 4445.4429 2776.2058
10 4178.2577 4445.4429
11 4928.2093 4178.2577
12 8432.8223 4928.2093
13 2727.7273 8432.8223
14 -325.8957 2727.7273
15 -3120.8494 -325.8957
16 -2084.1579 -3120.8494
17 1592.2167 -2084.1579
18 -381.4389 1592.2167
19 3709.3956 -381.4389
20 10849.8179 3709.3956
21 15923.6657 10849.8179
22 7385.2591 15923.6657
23 168.9894 7385.2591
24 7272.7499 168.9894
25 12366.9501 7272.7499
26 7807.3271 12366.9501
27 1681.6889 7807.3271
28 7503.4543 1681.6889
29 16883.5133 7503.4543
30 16413.0053 16883.5133
31 10903.1553 16413.0053
32 14544.9670 10903.1553
33 3341.1303 14544.9670
34 8666.1664 3341.1303
35 5542.6550 8666.1664
36 14326.8048 5542.6550
37 10386.3943 14326.8048
38 9657.5295 10386.3943
39 12543.6496 9657.5295
40 9294.3412 12543.6496
41 -3484.5057 9294.3412
42 -1901.0874 -3484.5057
43 -5897.9374 -1901.0874
44 -15438.2733 -5897.9374
45 -8950.8886 -15438.2733
46 -6125.4632 -8950.8886
47 -2920.5116 -6125.4632
48 -1151.5093 -2920.5116
49 -9014.0674 -1151.5093
50 -7404.3952 -9014.0674
51 -3987.8120 -7404.3952
52 -16750.2681 -3987.8120
53 -14920.0411 -16750.2681
54 -16745.1598 -14920.0411
55 -22780.0097 -16745.1598
56 -28673.8826 -22780.0097
57 -14759.3503 -28673.8826
58 -14104.2200 -14759.3503
59 -7719.3422 -14104.2200
60 -6749.5613 -7719.3422
61 -3123.5087 -6749.5613
62 1907.4052 -3123.5087
63 430.0261 1907.4052
64 5597.0316 430.0261
65 13476.2861 5597.0316
66 15363.7318 13476.2861
67 22525.2233 15363.7318
68 15941.1652 22525.2233
69 NA 15941.1652
> dum1 <- dum[2:length(myerror),]
> dum1
lag(myerror, k = 1) myerror
[1,] -13343.4955 -22131.3064
[2,] -11641.9709 -13343.4955
[3,] -7546.7033 -11641.9709
[4,] -3560.4011 -7546.7033
[5,] -13547.4692 -3560.4011
[6,] -12749.0510 -13547.4692
[7,] -8459.8272 -12749.0510
[8,] 2776.2058 -8459.8272
[9,] 4445.4429 2776.2058
[10,] 4178.2577 4445.4429
[11,] 4928.2093 4178.2577
[12,] 8432.8223 4928.2093
[13,] 2727.7273 8432.8223
[14,] -325.8957 2727.7273
[15,] -3120.8494 -325.8957
[16,] -2084.1579 -3120.8494
[17,] 1592.2167 -2084.1579
[18,] -381.4389 1592.2167
[19,] 3709.3956 -381.4389
[20,] 10849.8179 3709.3956
[21,] 15923.6657 10849.8179
[22,] 7385.2591 15923.6657
[23,] 168.9894 7385.2591
[24,] 7272.7499 168.9894
[25,] 12366.9501 7272.7499
[26,] 7807.3271 12366.9501
[27,] 1681.6889 7807.3271
[28,] 7503.4543 1681.6889
[29,] 16883.5133 7503.4543
[30,] 16413.0053 16883.5133
[31,] 10903.1553 16413.0053
[32,] 14544.9670 10903.1553
[33,] 3341.1303 14544.9670
[34,] 8666.1664 3341.1303
[35,] 5542.6550 8666.1664
[36,] 14326.8048 5542.6550
[37,] 10386.3943 14326.8048
[38,] 9657.5295 10386.3943
[39,] 12543.6496 9657.5295
[40,] 9294.3412 12543.6496
[41,] -3484.5057 9294.3412
[42,] -1901.0874 -3484.5057
[43,] -5897.9374 -1901.0874
[44,] -15438.2733 -5897.9374
[45,] -8950.8886 -15438.2733
[46,] -6125.4632 -8950.8886
[47,] -2920.5116 -6125.4632
[48,] -1151.5093 -2920.5116
[49,] -9014.0674 -1151.5093
[50,] -7404.3952 -9014.0674
[51,] -3987.8120 -7404.3952
[52,] -16750.2681 -3987.8120
[53,] -14920.0411 -16750.2681
[54,] -16745.1598 -14920.0411
[55,] -22780.0097 -16745.1598
[56,] -28673.8826 -22780.0097
[57,] -14759.3503 -28673.8826
[58,] -14104.2200 -14759.3503
[59,] -7719.3422 -14104.2200
[60,] -6749.5613 -7719.3422
[61,] -3123.5087 -6749.5613
[62,] 1907.4052 -3123.5087
[63,] 430.0261 1907.4052
[64,] 5597.0316 430.0261
[65,] 13476.2861 5597.0316
[66,] 15363.7318 13476.2861
[67,] 22525.2233 15363.7318
[68,] 15941.1652 22525.2233
> z <- as.data.frame(dum1)
> z
lag(myerror, k = 1) myerror
1 -13343.4955 -22131.3064
2 -11641.9709 -13343.4955
3 -7546.7033 -11641.9709
4 -3560.4011 -7546.7033
5 -13547.4692 -3560.4011
6 -12749.0510 -13547.4692
7 -8459.8272 -12749.0510
8 2776.2058 -8459.8272
9 4445.4429 2776.2058
10 4178.2577 4445.4429
11 4928.2093 4178.2577
12 8432.8223 4928.2093
13 2727.7273 8432.8223
14 -325.8957 2727.7273
15 -3120.8494 -325.8957
16 -2084.1579 -3120.8494
17 1592.2167 -2084.1579
18 -381.4389 1592.2167
19 3709.3956 -381.4389
20 10849.8179 3709.3956
21 15923.6657 10849.8179
22 7385.2591 15923.6657
23 168.9894 7385.2591
24 7272.7499 168.9894
25 12366.9501 7272.7499
26 7807.3271 12366.9501
27 1681.6889 7807.3271
28 7503.4543 1681.6889
29 16883.5133 7503.4543
30 16413.0053 16883.5133
31 10903.1553 16413.0053
32 14544.9670 10903.1553
33 3341.1303 14544.9670
34 8666.1664 3341.1303
35 5542.6550 8666.1664
36 14326.8048 5542.6550
37 10386.3943 14326.8048
38 9657.5295 10386.3943
39 12543.6496 9657.5295
40 9294.3412 12543.6496
41 -3484.5057 9294.3412
42 -1901.0874 -3484.5057
43 -5897.9374 -1901.0874
44 -15438.2733 -5897.9374
45 -8950.8886 -15438.2733
46 -6125.4632 -8950.8886
47 -2920.5116 -6125.4632
48 -1151.5093 -2920.5116
49 -9014.0674 -1151.5093
50 -7404.3952 -9014.0674
51 -3987.8120 -7404.3952
52 -16750.2681 -3987.8120
53 -14920.0411 -16750.2681
54 -16745.1598 -14920.0411
55 -22780.0097 -16745.1598
56 -28673.8826 -22780.0097
57 -14759.3503 -28673.8826
58 -14104.2200 -14759.3503
59 -7719.3422 -14104.2200
60 -6749.5613 -7719.3422
61 -3123.5087 -6749.5613
62 1907.4052 -3123.5087
63 430.0261 1907.4052
64 5597.0316 430.0261
65 13476.2861 5597.0316
66 15363.7318 13476.2861
67 22525.2233 15363.7318
68 15941.1652 22525.2233
> 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/7p9v11258715003.ps",horizontal=F,pagecentre=F,paper="special",width=8.3333333333333,height=5.5555555555556)
> acf(mysum$resid, lag.max=length(mysum$resid)/2, main='Residual Autocorrelation Function')
> grid()
> dev.off()
null device
1
> postscript(file="/var/www/html/rcomp/tmp/8pjxg1258715003.ps",horizontal=F,pagecentre=F,paper="special",width=8.3333333333333,height=5.5555555555556)
> pacf(mysum$resid, lag.max=length(mysum$resid)/2, main='Residual Partial Autocorrelation Function')
> grid()
> dev.off()
null device
1
> postscript(file="/var/www/html/rcomp/tmp/9h8mq1258715003.ps",horizontal=F,pagecentre=F,paper="special",width=8.3333333333333,height=5.5555555555556)
> opar <- par(mfrow = c(2,2), oma = c(0, 0, 1.1, 0))
> plot(mylm, las = 1, sub='Residual Diagnostics')
> par(opar)
> dev.off()
null device
1
> if (n > n25) {
+ postscript(file="/var/www/html/rcomp/tmp/10469u1258715003.ps",horizontal=F,pagecentre=F,paper="special",width=8.3333333333333,height=5.5555555555556)
+ plot(kp3:nmkm3,gqarr[,2], main='Goldfeld-Quandt test',ylab='2-sided p-value',xlab='breakpoint')
+ grid()
+ dev.off()
+ }
null device
1
>
> #Note: the /var/www/html/rcomp/createtable file can be downloaded at http://www.wessa.net/cretab
> load(file="/var/www/html/rcomp/createtable")
>
> a<-table.start()
> a<-table.row.start(a)
> a<-table.element(a, 'Multiple Linear Regression - Estimated Regression Equation', 1, TRUE)
> a<-table.row.end(a)
> myeq <- colnames(x)[1]
> myeq <- paste(myeq, '[t] = ', sep='')
> for (i in 1:k){
+ if (mysum$coefficients[i,1] > 0) myeq <- paste(myeq, '+', '')
+ myeq <- paste(myeq, mysum$coefficients[i,1], sep=' ')
+ if (rownames(mysum$coefficients)[i] != '(Intercept)') {
+ myeq <- paste(myeq, rownames(mysum$coefficients)[i], sep='')
+ if (rownames(mysum$coefficients)[i] != 't') myeq <- paste(myeq, '[t]', sep='')
+ }
+ }
> myeq <- paste(myeq, ' + e[t]')
> a<-table.row.start(a)
> a<-table.element(a, myeq)
> a<-table.row.end(a)
> a<-table.end(a)
> table.save(a,file="/var/www/html/rcomp/tmp/11s5181258715003.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/12lpp91258715003.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/13iskh1258715003.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/148lin1258715003.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/15673a1258715003.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/1692jm1258715003.tab")
+ }
>
> system("convert tmp/13vcg1258715003.ps tmp/13vcg1258715003.png")
> system("convert tmp/2iehl1258715003.ps tmp/2iehl1258715003.png")
> system("convert tmp/3leyy1258715003.ps tmp/3leyy1258715003.png")
> system("convert tmp/4046s1258715003.ps tmp/4046s1258715003.png")
> system("convert tmp/5f02p1258715003.ps tmp/5f02p1258715003.png")
> system("convert tmp/6j5b91258715003.ps tmp/6j5b91258715003.png")
> system("convert tmp/7p9v11258715003.ps tmp/7p9v11258715003.png")
> system("convert tmp/8pjxg1258715003.ps tmp/8pjxg1258715003.png")
> system("convert tmp/9h8mq1258715003.ps tmp/9h8mq1258715003.png")
> system("convert tmp/10469u1258715003.ps tmp/10469u1258715003.png")
>
>
> proc.time()
user system elapsed
2.500 1.551 3.016