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(577992,0,565464,0,547344,0,554788,0,562325,0,560854,0,555332,0,543599,0,536662,0,542722,0,593530,1,610763,1,612613,1,611324,1,594167,1,595454,1,590865,1,589379,1,584428,1,573100,1,567456,1,569028,1,620735,1,628884,1,628232,1,612117,1,595404,1,597141,1,593408,1,590072,1,579799,1,574205,1,572775,1,572942,1,619567,1,625809,1,619916,1,587625,1,565742,1,557274,1,560576,1,548854,1,531673,1,525919,1,511038,1,498662,1,555362,1,564591,1,541657,1,527070,1,509846,1,514258,1,516922,1,507561,1,492622,1,490243,1,469357,1,477580,1,528379,1,533590,1,517945,1,506174,1,501866,1,516141,1,528222,1,532638,1,536322,1,536535,1,523597,1,536214,1,586570,1,596594,1,580523,1),dim=c(2,73),dimnames=list(c('Y','X'),1:73))
> y <- array(NA,dim=c(2,73),dimnames=list(c('Y','X'),1:73))
> 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 577992 0 1 0 0 0 0 0 0 0 0 0 0 1
2 565464 0 0 1 0 0 0 0 0 0 0 0 0 2
3 547344 0 0 0 1 0 0 0 0 0 0 0 0 3
4 554788 0 0 0 0 1 0 0 0 0 0 0 0 4
5 562325 0 0 0 0 0 1 0 0 0 0 0 0 5
6 560854 0 0 0 0 0 0 1 0 0 0 0 0 6
7 555332 0 0 0 0 0 0 0 1 0 0 0 0 7
8 543599 0 0 0 0 0 0 0 0 1 0 0 0 8
9 536662 0 0 0 0 0 0 0 0 0 1 0 0 9
10 542722 0 0 0 0 0 0 0 0 0 0 1 0 10
11 593530 1 0 0 0 0 0 0 0 0 0 0 1 11
12 610763 1 0 0 0 0 0 0 0 0 0 0 0 12
13 612613 1 1 0 0 0 0 0 0 0 0 0 0 13
14 611324 1 0 1 0 0 0 0 0 0 0 0 0 14
15 594167 1 0 0 1 0 0 0 0 0 0 0 0 15
16 595454 1 0 0 0 1 0 0 0 0 0 0 0 16
17 590865 1 0 0 0 0 1 0 0 0 0 0 0 17
18 589379 1 0 0 0 0 0 1 0 0 0 0 0 18
19 584428 1 0 0 0 0 0 0 1 0 0 0 0 19
20 573100 1 0 0 0 0 0 0 0 1 0 0 0 20
21 567456 1 0 0 0 0 0 0 0 0 1 0 0 21
22 569028 1 0 0 0 0 0 0 0 0 0 1 0 22
23 620735 1 0 0 0 0 0 0 0 0 0 0 1 23
24 628884 1 0 0 0 0 0 0 0 0 0 0 0 24
25 628232 1 1 0 0 0 0 0 0 0 0 0 0 25
26 612117 1 0 1 0 0 0 0 0 0 0 0 0 26
27 595404 1 0 0 1 0 0 0 0 0 0 0 0 27
28 597141 1 0 0 0 1 0 0 0 0 0 0 0 28
29 593408 1 0 0 0 0 1 0 0 0 0 0 0 29
30 590072 1 0 0 0 0 0 1 0 0 0 0 0 30
31 579799 1 0 0 0 0 0 0 1 0 0 0 0 31
32 574205 1 0 0 0 0 0 0 0 1 0 0 0 32
33 572775 1 0 0 0 0 0 0 0 0 1 0 0 33
34 572942 1 0 0 0 0 0 0 0 0 0 1 0 34
35 619567 1 0 0 0 0 0 0 0 0 0 0 1 35
36 625809 1 0 0 0 0 0 0 0 0 0 0 0 36
37 619916 1 1 0 0 0 0 0 0 0 0 0 0 37
38 587625 1 0 1 0 0 0 0 0 0 0 0 0 38
39 565742 1 0 0 1 0 0 0 0 0 0 0 0 39
40 557274 1 0 0 0 1 0 0 0 0 0 0 0 40
41 560576 1 0 0 0 0 1 0 0 0 0 0 0 41
42 548854 1 0 0 0 0 0 1 0 0 0 0 0 42
43 531673 1 0 0 0 0 0 0 1 0 0 0 0 43
44 525919 1 0 0 0 0 0 0 0 1 0 0 0 44
45 511038 1 0 0 0 0 0 0 0 0 1 0 0 45
46 498662 1 0 0 0 0 0 0 0 0 0 1 0 46
47 555362 1 0 0 0 0 0 0 0 0 0 0 1 47
48 564591 1 0 0 0 0 0 0 0 0 0 0 0 48
49 541657 1 1 0 0 0 0 0 0 0 0 0 0 49
50 527070 1 0 1 0 0 0 0 0 0 0 0 0 50
51 509846 1 0 0 1 0 0 0 0 0 0 0 0 51
52 514258 1 0 0 0 1 0 0 0 0 0 0 0 52
53 516922 1 0 0 0 0 1 0 0 0 0 0 0 53
54 507561 1 0 0 0 0 0 1 0 0 0 0 0 54
55 492622 1 0 0 0 0 0 0 1 0 0 0 0 55
56 490243 1 0 0 0 0 0 0 0 1 0 0 0 56
57 469357 1 0 0 0 0 0 0 0 0 1 0 0 57
58 477580 1 0 0 0 0 0 0 0 0 0 1 0 58
59 528379 1 0 0 0 0 0 0 0 0 0 0 1 59
60 533590 1 0 0 0 0 0 0 0 0 0 0 0 60
61 517945 1 1 0 0 0 0 0 0 0 0 0 0 61
62 506174 1 0 1 0 0 0 0 0 0 0 0 0 62
63 501866 1 0 0 1 0 0 0 0 0 0 0 0 63
64 516141 1 0 0 0 1 0 0 0 0 0 0 0 64
65 528222 1 0 0 0 0 1 0 0 0 0 0 0 65
66 532638 1 0 0 0 0 0 1 0 0 0 0 0 66
67 536322 1 0 0 0 0 0 0 1 0 0 0 0 67
68 536535 1 0 0 0 0 0 0 0 1 0 0 0 68
69 523597 1 0 0 0 0 0 0 0 0 1 0 0 69
70 536214 1 0 0 0 0 0 0 0 0 0 1 0 70
71 586570 1 0 0 0 0 0 0 0 0 0 0 1 71
72 596594 1 0 0 0 0 0 0 0 0 0 0 0 72
73 580523 1 1 0 0 0 0 0 0 0 0 0 0 73
> k <- length(x[1,])
> df <- as.data.frame(x)
> (mylm <- lm(df))
Call:
lm(formula = df)
Coefficients:
(Intercept) X M1 M2 M3 M4
604336 52211 -10737 -31416 -45813 -40861
M5 M6 M7 M8 M9 M10
-36479 -38802 -45495 -50086 -59035 -54820
M11 t
-10852 -1504
> (mysum <- summary(mylm))
Call:
lm(formula = df)
Residuals:
Min 1Q Median 3Q Max
-42417 -19637 1754 16961 48347
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 604335.5 13504.1 44.752 < 2e-16 ***
X 52211.4 10825.4 4.823 1.03e-05 ***
M1 -10737.1 14061.2 -0.764 0.448152
M2 -31416.0 14623.0 -2.148 0.035799 *
M3 -45812.6 14615.9 -3.134 0.002682 **
M4 -40860.6 14610.8 -2.797 0.006961 **
M5 -36479.5 14607.9 -2.497 0.015326 *
M6 -38801.9 14607.0 -2.656 0.010143 *
M7 -45494.8 14608.3 -3.114 0.002844 **
M8 -50086.4 14611.7 -3.428 0.001115 **
M9 -59034.9 14617.2 -4.039 0.000157 ***
M10 -54820.3 14624.8 -3.748 0.000407 ***
M11 -10852.2 14535.8 -0.747 0.458280
t -1504.2 175.6 -8.565 6.11e-12 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 25170 on 59 degrees of freedom
Multiple R-squared: 0.6652, Adjusted R-squared: 0.5915
F-statistic: 9.019 on 13 and 59 DF, p-value: 8.287e-10
> 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,] 1.494797e-02 2.989595e-02 0.985052026
[2,] 5.136936e-03 1.027387e-02 0.994863064
[3,] 1.443920e-03 2.887840e-03 0.998556080
[4,] 3.584122e-04 7.168244e-04 0.999641588
[5,] 7.467641e-05 1.493528e-04 0.999925324
[6,] 2.167145e-05 4.334291e-05 0.999978329
[7,] 3.698906e-06 7.397812e-06 0.999996301
[8,] 8.560054e-07 1.712011e-06 0.999999144
[9,] 1.891737e-07 3.783474e-07 0.999999811
[10,] 1.338179e-07 2.676358e-07 0.999999866
[11,] 4.489302e-08 8.978603e-08 0.999999955
[12,] 1.733244e-08 3.466488e-08 0.999999983
[13,] 1.139722e-08 2.279444e-08 0.999999989
[14,] 7.253080e-09 1.450616e-08 0.999999993
[15,] 8.114547e-09 1.622909e-08 0.999999992
[16,] 2.751932e-09 5.503864e-09 0.999999997
[17,] 8.623498e-10 1.724700e-09 0.999999999
[18,] 3.452089e-10 6.904179e-10 1.000000000
[19,] 1.176574e-10 2.353149e-10 1.000000000
[20,] 4.653287e-11 9.306574e-11 1.000000000
[21,] 5.777776e-11 1.155555e-10 1.000000000
[22,] 9.477227e-09 1.895445e-08 0.999999991
[23,] 4.065177e-07 8.130355e-07 0.999999593
[24,] 1.723594e-05 3.447189e-05 0.999982764
[25,] 1.122835e-04 2.245671e-04 0.999887716
[26,] 8.260270e-04 1.652054e-03 0.999173973
[27,] 4.819737e-03 9.639475e-03 0.995180263
[28,] 1.127795e-02 2.255589e-02 0.988722054
[29,] 3.674451e-02 7.348901e-02 0.963255494
[30,] 8.299258e-02 1.659852e-01 0.917007416
[31,] 9.383191e-02 1.876638e-01 0.906168092
[32,] 1.261167e-01 2.522335e-01 0.873883257
[33,] 2.071987e-01 4.143975e-01 0.792801272
[34,] 4.302225e-01 8.604449e-01 0.569777540
[35,] 6.491632e-01 7.016736e-01 0.350836810
[36,] 8.298823e-01 3.402355e-01 0.170117731
[37,] 9.523938e-01 9.521234e-02 0.047606171
[38,] 9.925769e-01 1.484628e-02 0.007423142
[39,] 9.919865e-01 1.602697e-02 0.008013485
[40,] 9.951238e-01 9.752349e-03 0.004876175
> postscript(file="/var/www/html/rcomp/tmp/10yt21260651763.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/24nq51260651763.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/3nph11260651763.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/4rmmw1260651763.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/571si1260651763.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 = 73
Frequency = 1
1 2 3 4 5 6
-14102.3235 -4447.2418 -6666.4085 -2670.2418 1989.7582 4345.4248
7 8 9 10 11 12
7020.4248 1383.2582 4898.9248 8248.4248 -35618.9358 -27733.9358
13 14 15 16 17 18
-13642.7152 7251.3664 5995.1998 3834.3664 -3631.6336 -1290.9669
19 20 21 22 23 24
1955.0331 -3277.1336 1531.5331 393.0331 9636.1052 8437.1052
25 26 27 28 29 30
20026.3258 26094.4074 25282.2407 23571.4074 16961.4074 17452.0741
31 32 33 34 35 36
15376.0741 15877.9074 24900.5741 22357.0741 26518.1462 23412.1462
37 38 39 40 41 42
29760.3668 19652.4484 13670.2817 1754.4484 2179.4484 -5715.8850
43 44 45 46 47 48
-14699.8850 -14358.0516 -18786.3850 -33872.8850 -19636.8128 -19755.8128
49 50 51 52 53 54
-30448.5923 -22852.5107 -24175.6773 -23211.5107 -23424.5107 -28958.8440
55 56 57 58 59 60
-35700.8440 -31984.0107 -42417.3440 -36904.8440 -28569.7719 -32706.7719
61 62 63 64 65 66
-36110.5513 -25698.4697 -14105.6364 -3278.4697 5925.5303 14168.1970
67 68 69 70 71 72
26049.1970 32358.0303 29872.6970 39779.1970 47671.2691 48347.2691
73
44517.4897
> postscript(file="/var/www/html/rcomp/tmp/6na4v1260651763.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 = 73
Frequency = 1
lag(myerror, k = 1) myerror
0 -14102.3235 NA
1 -4447.2418 -14102.3235
2 -6666.4085 -4447.2418
3 -2670.2418 -6666.4085
4 1989.7582 -2670.2418
5 4345.4248 1989.7582
6 7020.4248 4345.4248
7 1383.2582 7020.4248
8 4898.9248 1383.2582
9 8248.4248 4898.9248
10 -35618.9358 8248.4248
11 -27733.9358 -35618.9358
12 -13642.7152 -27733.9358
13 7251.3664 -13642.7152
14 5995.1998 7251.3664
15 3834.3664 5995.1998
16 -3631.6336 3834.3664
17 -1290.9669 -3631.6336
18 1955.0331 -1290.9669
19 -3277.1336 1955.0331
20 1531.5331 -3277.1336
21 393.0331 1531.5331
22 9636.1052 393.0331
23 8437.1052 9636.1052
24 20026.3258 8437.1052
25 26094.4074 20026.3258
26 25282.2407 26094.4074
27 23571.4074 25282.2407
28 16961.4074 23571.4074
29 17452.0741 16961.4074
30 15376.0741 17452.0741
31 15877.9074 15376.0741
32 24900.5741 15877.9074
33 22357.0741 24900.5741
34 26518.1462 22357.0741
35 23412.1462 26518.1462
36 29760.3668 23412.1462
37 19652.4484 29760.3668
38 13670.2817 19652.4484
39 1754.4484 13670.2817
40 2179.4484 1754.4484
41 -5715.8850 2179.4484
42 -14699.8850 -5715.8850
43 -14358.0516 -14699.8850
44 -18786.3850 -14358.0516
45 -33872.8850 -18786.3850
46 -19636.8128 -33872.8850
47 -19755.8128 -19636.8128
48 -30448.5923 -19755.8128
49 -22852.5107 -30448.5923
50 -24175.6773 -22852.5107
51 -23211.5107 -24175.6773
52 -23424.5107 -23211.5107
53 -28958.8440 -23424.5107
54 -35700.8440 -28958.8440
55 -31984.0107 -35700.8440
56 -42417.3440 -31984.0107
57 -36904.8440 -42417.3440
58 -28569.7719 -36904.8440
59 -32706.7719 -28569.7719
60 -36110.5513 -32706.7719
61 -25698.4697 -36110.5513
62 -14105.6364 -25698.4697
63 -3278.4697 -14105.6364
64 5925.5303 -3278.4697
65 14168.1970 5925.5303
66 26049.1970 14168.1970
67 32358.0303 26049.1970
68 29872.6970 32358.0303
69 39779.1970 29872.6970
70 47671.2691 39779.1970
71 48347.2691 47671.2691
72 44517.4897 48347.2691
73 NA 44517.4897
> dum1 <- dum[2:length(myerror),]
> dum1
lag(myerror, k = 1) myerror
[1,] -4447.2418 -14102.3235
[2,] -6666.4085 -4447.2418
[3,] -2670.2418 -6666.4085
[4,] 1989.7582 -2670.2418
[5,] 4345.4248 1989.7582
[6,] 7020.4248 4345.4248
[7,] 1383.2582 7020.4248
[8,] 4898.9248 1383.2582
[9,] 8248.4248 4898.9248
[10,] -35618.9358 8248.4248
[11,] -27733.9358 -35618.9358
[12,] -13642.7152 -27733.9358
[13,] 7251.3664 -13642.7152
[14,] 5995.1998 7251.3664
[15,] 3834.3664 5995.1998
[16,] -3631.6336 3834.3664
[17,] -1290.9669 -3631.6336
[18,] 1955.0331 -1290.9669
[19,] -3277.1336 1955.0331
[20,] 1531.5331 -3277.1336
[21,] 393.0331 1531.5331
[22,] 9636.1052 393.0331
[23,] 8437.1052 9636.1052
[24,] 20026.3258 8437.1052
[25,] 26094.4074 20026.3258
[26,] 25282.2407 26094.4074
[27,] 23571.4074 25282.2407
[28,] 16961.4074 23571.4074
[29,] 17452.0741 16961.4074
[30,] 15376.0741 17452.0741
[31,] 15877.9074 15376.0741
[32,] 24900.5741 15877.9074
[33,] 22357.0741 24900.5741
[34,] 26518.1462 22357.0741
[35,] 23412.1462 26518.1462
[36,] 29760.3668 23412.1462
[37,] 19652.4484 29760.3668
[38,] 13670.2817 19652.4484
[39,] 1754.4484 13670.2817
[40,] 2179.4484 1754.4484
[41,] -5715.8850 2179.4484
[42,] -14699.8850 -5715.8850
[43,] -14358.0516 -14699.8850
[44,] -18786.3850 -14358.0516
[45,] -33872.8850 -18786.3850
[46,] -19636.8128 -33872.8850
[47,] -19755.8128 -19636.8128
[48,] -30448.5923 -19755.8128
[49,] -22852.5107 -30448.5923
[50,] -24175.6773 -22852.5107
[51,] -23211.5107 -24175.6773
[52,] -23424.5107 -23211.5107
[53,] -28958.8440 -23424.5107
[54,] -35700.8440 -28958.8440
[55,] -31984.0107 -35700.8440
[56,] -42417.3440 -31984.0107
[57,] -36904.8440 -42417.3440
[58,] -28569.7719 -36904.8440
[59,] -32706.7719 -28569.7719
[60,] -36110.5513 -32706.7719
[61,] -25698.4697 -36110.5513
[62,] -14105.6364 -25698.4697
[63,] -3278.4697 -14105.6364
[64,] 5925.5303 -3278.4697
[65,] 14168.1970 5925.5303
[66,] 26049.1970 14168.1970
[67,] 32358.0303 26049.1970
[68,] 29872.6970 32358.0303
[69,] 39779.1970 29872.6970
[70,] 47671.2691 39779.1970
[71,] 48347.2691 47671.2691
[72,] 44517.4897 48347.2691
> z <- as.data.frame(dum1)
> z
lag(myerror, k = 1) myerror
1 -4447.2418 -14102.3235
2 -6666.4085 -4447.2418
3 -2670.2418 -6666.4085
4 1989.7582 -2670.2418
5 4345.4248 1989.7582
6 7020.4248 4345.4248
7 1383.2582 7020.4248
8 4898.9248 1383.2582
9 8248.4248 4898.9248
10 -35618.9358 8248.4248
11 -27733.9358 -35618.9358
12 -13642.7152 -27733.9358
13 7251.3664 -13642.7152
14 5995.1998 7251.3664
15 3834.3664 5995.1998
16 -3631.6336 3834.3664
17 -1290.9669 -3631.6336
18 1955.0331 -1290.9669
19 -3277.1336 1955.0331
20 1531.5331 -3277.1336
21 393.0331 1531.5331
22 9636.1052 393.0331
23 8437.1052 9636.1052
24 20026.3258 8437.1052
25 26094.4074 20026.3258
26 25282.2407 26094.4074
27 23571.4074 25282.2407
28 16961.4074 23571.4074
29 17452.0741 16961.4074
30 15376.0741 17452.0741
31 15877.9074 15376.0741
32 24900.5741 15877.9074
33 22357.0741 24900.5741
34 26518.1462 22357.0741
35 23412.1462 26518.1462
36 29760.3668 23412.1462
37 19652.4484 29760.3668
38 13670.2817 19652.4484
39 1754.4484 13670.2817
40 2179.4484 1754.4484
41 -5715.8850 2179.4484
42 -14699.8850 -5715.8850
43 -14358.0516 -14699.8850
44 -18786.3850 -14358.0516
45 -33872.8850 -18786.3850
46 -19636.8128 -33872.8850
47 -19755.8128 -19636.8128
48 -30448.5923 -19755.8128
49 -22852.5107 -30448.5923
50 -24175.6773 -22852.5107
51 -23211.5107 -24175.6773
52 -23424.5107 -23211.5107
53 -28958.8440 -23424.5107
54 -35700.8440 -28958.8440
55 -31984.0107 -35700.8440
56 -42417.3440 -31984.0107
57 -36904.8440 -42417.3440
58 -28569.7719 -36904.8440
59 -32706.7719 -28569.7719
60 -36110.5513 -32706.7719
61 -25698.4697 -36110.5513
62 -14105.6364 -25698.4697
63 -3278.4697 -14105.6364
64 5925.5303 -3278.4697
65 14168.1970 5925.5303
66 26049.1970 14168.1970
67 32358.0303 26049.1970
68 29872.6970 32358.0303
69 39779.1970 29872.6970
70 47671.2691 39779.1970
71 48347.2691 47671.2691
72 44517.4897 48347.2691
> 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/7h67e1260651763.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/8xe4s1260651763.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/96l8r1260651763.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/10i5z71260651763.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/11bnvb1260651763.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/12zz9o1260651763.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/13rw8x1260651763.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/14jl2y1260651764.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/152zin1260651764.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/16gq5a1260651764.tab")
+ }
> try(system("convert tmp/10yt21260651763.ps tmp/10yt21260651763.png",intern=TRUE))
character(0)
> try(system("convert tmp/24nq51260651763.ps tmp/24nq51260651763.png",intern=TRUE))
character(0)
> try(system("convert tmp/3nph11260651763.ps tmp/3nph11260651763.png",intern=TRUE))
character(0)
> try(system("convert tmp/4rmmw1260651763.ps tmp/4rmmw1260651763.png",intern=TRUE))
character(0)
> try(system("convert tmp/571si1260651763.ps tmp/571si1260651763.png",intern=TRUE))
character(0)
> try(system("convert tmp/6na4v1260651763.ps tmp/6na4v1260651763.png",intern=TRUE))
character(0)
> try(system("convert tmp/7h67e1260651763.ps tmp/7h67e1260651763.png",intern=TRUE))
character(0)
> try(system("convert tmp/8xe4s1260651763.ps tmp/8xe4s1260651763.png",intern=TRUE))
character(0)
> try(system("convert tmp/96l8r1260651763.ps tmp/96l8r1260651763.png",intern=TRUE))
character(0)
> try(system("convert tmp/10i5z71260651763.ps tmp/10i5z71260651763.png",intern=TRUE))
character(0)
>
>
> proc.time()
user system elapsed
2.559 1.545 4.161