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(2440.25
+ ,0
+ ,2350.44
+ ,2408.64
+ ,0
+ ,2440.25
+ ,2472.81
+ ,0
+ ,2408.64
+ ,2407.6
+ ,0
+ ,2472.81
+ ,2454.62
+ ,0
+ ,2407.6
+ ,2448.05
+ ,0
+ ,2454.62
+ ,2497.84
+ ,0
+ ,2448.05
+ ,2645.64
+ ,0
+ ,2497.84
+ ,2756.76
+ ,0
+ ,2645.64
+ ,2849.27
+ ,0
+ ,2756.76
+ ,2921.44
+ ,0
+ ,2849.27
+ ,2981.85
+ ,0
+ ,2921.44
+ ,3080.58
+ ,0
+ ,2981.85
+ ,3106.22
+ ,0
+ ,3080.58
+ ,3119.31
+ ,0
+ ,3106.22
+ ,3061.26
+ ,0
+ ,3119.31
+ ,3097.31
+ ,0
+ ,3061.26
+ ,3161.69
+ ,0
+ ,3097.31
+ ,3257.16
+ ,0
+ ,3161.69
+ ,3277.01
+ ,0
+ ,3257.16
+ ,3295.32
+ ,0
+ ,3277.01
+ ,3363.99
+ ,0
+ ,3295.32
+ ,3494.17
+ ,0
+ ,3363.99
+ ,3667.03
+ ,0
+ ,3494.17
+ ,3813.06
+ ,0
+ ,3667.03
+ ,3917.96
+ ,0
+ ,3813.06
+ ,3895.51
+ ,0
+ ,3917.96
+ ,3801.06
+ ,0
+ ,3895.51
+ ,3570.12
+ ,0
+ ,3801.06
+ ,3701.61
+ ,0
+ ,3570.12
+ ,3862.27
+ ,0
+ ,3701.61
+ ,3970.1
+ ,0
+ ,3862.27
+ ,4138.52
+ ,0
+ ,3970.1
+ ,4199.75
+ ,0
+ ,4138.52
+ ,4290.89
+ ,0
+ ,4199.75
+ ,4443.91
+ ,0
+ ,4290.89
+ ,4502.64
+ ,0
+ ,4443.91
+ ,4356.98
+ ,0
+ ,4502.64
+ ,4591.27
+ ,0
+ ,4356.98
+ ,4696.96
+ ,0
+ ,4591.27
+ ,4621.4
+ ,0
+ ,4696.96
+ ,4562.84
+ ,0
+ ,4621.4
+ ,4202.52
+ ,0
+ ,4562.84
+ ,4296.49
+ ,0
+ ,4202.52
+ ,4435.23
+ ,0
+ ,4296.49
+ ,4105.18
+ ,0
+ ,4435.23
+ ,4116.68
+ ,0
+ ,4105.18
+ ,3844.49
+ ,0
+ ,4116.68
+ ,3720.98
+ ,0
+ ,3844.49
+ ,3674.4
+ ,0
+ ,3720.98
+ ,3857.62
+ ,0
+ ,3674.4
+ ,3801.06
+ ,0
+ ,3857.62
+ ,3504.37
+ ,0
+ ,3801.06
+ ,3032.6
+ ,0
+ ,3504.37
+ ,3047.03
+ ,0
+ ,3032.6
+ ,2962.34
+ ,1
+ ,3047.03
+ ,2197.82
+ ,1
+ ,2962.34
+ ,2014.45
+ ,1
+ ,2197.82
+ ,1862.83
+ ,1
+ ,2014.45
+ ,1905.41
+ ,1
+ ,1862.83
+ ,1810.99
+ ,1
+ ,1905.41
+ ,1670.07
+ ,1
+ ,1810.99
+ ,1864.44
+ ,1
+ ,1670.07
+ ,2052.02
+ ,1
+ ,1864.44
+ ,2029.6
+ ,1
+ ,2052.02
+ ,2070.83
+ ,1
+ ,2029.6
+ ,2293.41
+ ,1
+ ,2070.83
+ ,2443.27
+ ,1
+ ,2293.41
+ ,2513.17
+ ,1
+ ,2443.27
+ ,2466.92
+ ,1
+ ,2513.17)
+ ,dim=c(3
+ ,70)
+ ,dimnames=list(c('Y'
+ ,'X'
+ ,'Y1')
+ ,1:70))
> y <- array(NA,dim=c(3,70),dimnames=list(c('Y','X','Y1'),1:70))
> 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 = 'Do not include Seasonal 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 Y1 t
1 2440.25 0 2350.44 1
2 2408.64 0 2440.25 2
3 2472.81 0 2408.64 3
4 2407.60 0 2472.81 4
5 2454.62 0 2407.60 5
6 2448.05 0 2454.62 6
7 2497.84 0 2448.05 7
8 2645.64 0 2497.84 8
9 2756.76 0 2645.64 9
10 2849.27 0 2756.76 10
11 2921.44 0 2849.27 11
12 2981.85 0 2921.44 12
13 3080.58 0 2981.85 13
14 3106.22 0 3080.58 14
15 3119.31 0 3106.22 15
16 3061.26 0 3119.31 16
17 3097.31 0 3061.26 17
18 3161.69 0 3097.31 18
19 3257.16 0 3161.69 19
20 3277.01 0 3257.16 20
21 3295.32 0 3277.01 21
22 3363.99 0 3295.32 22
23 3494.17 0 3363.99 23
24 3667.03 0 3494.17 24
25 3813.06 0 3667.03 25
26 3917.96 0 3813.06 26
27 3895.51 0 3917.96 27
28 3801.06 0 3895.51 28
29 3570.12 0 3801.06 29
30 3701.61 0 3570.12 30
31 3862.27 0 3701.61 31
32 3970.10 0 3862.27 32
33 4138.52 0 3970.10 33
34 4199.75 0 4138.52 34
35 4290.89 0 4199.75 35
36 4443.91 0 4290.89 36
37 4502.64 0 4443.91 37
38 4356.98 0 4502.64 38
39 4591.27 0 4356.98 39
40 4696.96 0 4591.27 40
41 4621.40 0 4696.96 41
42 4562.84 0 4621.40 42
43 4202.52 0 4562.84 43
44 4296.49 0 4202.52 44
45 4435.23 0 4296.49 45
46 4105.18 0 4435.23 46
47 4116.68 0 4105.18 47
48 3844.49 0 4116.68 48
49 3720.98 0 3844.49 49
50 3674.40 0 3720.98 50
51 3857.62 0 3674.40 51
52 3801.06 0 3857.62 52
53 3504.37 0 3801.06 53
54 3032.60 0 3504.37 54
55 3047.03 0 3032.60 55
56 2962.34 1 3047.03 56
57 2197.82 1 2962.34 57
58 2014.45 1 2197.82 58
59 1862.83 1 2014.45 59
60 1905.41 1 1862.83 60
61 1810.99 1 1905.41 61
62 1670.07 1 1810.99 62
63 1864.44 1 1670.07 63
64 2052.02 1 1864.44 64
65 2029.60 1 2052.02 65
66 2070.83 1 2029.60 66
67 2293.41 1 2070.83 67
68 2443.27 1 2293.41 68
69 2513.17 1 2443.27 69
70 2466.92 1 2513.17 70
> k <- length(x[1,])
> df <- as.data.frame(x)
> (mylm <- lm(df))
Call:
lm(formula = df)
Coefficients:
(Intercept) X Y1 t
227.0398 -97.2525 0.9455 -0.7861
> (mysum <- summary(mylm))
Call:
lm(formula = df)
Residuals:
Min 1Q Median 3Q Max
-688.066 -65.536 5.529 113.839 275.344
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 227.0399 132.0454 1.719 0.0902 .
X -97.2525 139.7743 -0.696 0.4890
Y1 0.9455 0.0480 19.698 <2e-16 ***
t -0.7861 2.1460 -0.366 0.7153
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 167.1 on 66 degrees of freedom
Multiple R-squared: 0.9628, Adjusted R-squared: 0.9611
F-statistic: 568.7 on 3 and 66 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,] 4.369972e-03 8.739944e-03 0.9956300
[2,] 4.470491e-02 8.940982e-02 0.9552951
[3,] 3.582796e-02 7.165593e-02 0.9641720
[4,] 1.383104e-02 2.766208e-02 0.9861690
[5,] 4.843815e-03 9.687630e-03 0.9951562
[6,] 1.626840e-03 3.253679e-03 0.9983732
[7,] 5.191990e-04 1.038398e-03 0.9994808
[8,] 2.096952e-04 4.193904e-04 0.9997903
[9,] 9.198986e-05 1.839797e-04 0.9999080
[10,] 1.534390e-04 3.068781e-04 0.9998466
[11,] 5.642143e-05 1.128429e-04 0.9999436
[12,] 1.774805e-05 3.549610e-05 0.9999823
[13,] 5.867288e-06 1.173458e-05 0.9999941
[14,] 2.156874e-06 4.313748e-06 0.9999978
[15,] 8.136162e-07 1.627232e-06 0.9999992
[16,] 2.342808e-07 4.685617e-07 0.9999998
[17,] 1.201232e-07 2.402463e-07 0.9999999
[18,] 1.786012e-07 3.572024e-07 0.9999998
[19,] 1.239152e-07 2.478305e-07 0.9999999
[20,] 4.066552e-08 8.133105e-08 1.0000000
[21,] 2.697976e-08 5.395951e-08 1.0000000
[22,] 7.934359e-08 1.586872e-07 0.9999999
[23,] 1.550584e-05 3.101167e-05 0.9999845
[24,] 6.805784e-06 1.361157e-05 0.9999932
[25,] 3.463860e-06 6.927720e-06 0.9999965
[26,] 1.431727e-06 2.863453e-06 0.9999986
[27,] 1.017852e-06 2.035704e-06 0.9999990
[28,] 3.868120e-07 7.736240e-07 0.9999996
[29,] 1.652942e-07 3.305884e-07 0.9999998
[30,] 1.499335e-07 2.998670e-07 0.9999999
[31,] 6.614125e-08 1.322825e-07 0.9999999
[32,] 1.155835e-07 2.311670e-07 0.9999999
[33,] 4.732078e-07 9.464157e-07 0.9999995
[34,] 5.821634e-07 1.164327e-06 0.9999994
[35,] 4.821131e-07 9.642262e-07 0.9999995
[36,] 5.250446e-07 1.050089e-06 0.9999995
[37,] 6.838065e-05 1.367613e-04 0.9999316
[38,] 8.965965e-05 1.793193e-04 0.9999103
[39,] 3.269691e-04 6.539381e-04 0.9996730
[40,] 2.525709e-03 5.051418e-03 0.9974743
[41,] 3.692723e-03 7.385446e-03 0.9963073
[42,] 5.565967e-03 1.113193e-02 0.9944340
[43,] 3.795075e-03 7.590150e-03 0.9962049
[44,] 2.860568e-03 5.721137e-03 0.9971394
[45,] 2.210554e-02 4.421108e-02 0.9778945
[46,] 4.391943e-02 8.783886e-02 0.9560806
[47,] 4.915525e-02 9.831050e-02 0.9508447
[48,] 1.121181e-01 2.242362e-01 0.8878819
[49,] 8.845594e-02 1.769119e-01 0.9115441
[50,] 8.268760e-01 3.462480e-01 0.1731240
[51,] 8.771669e-01 2.456662e-01 0.1228331
[52,] 8.530488e-01 2.939024e-01 0.1469512
[53,] 7.919925e-01 4.160151e-01 0.2080075
[54,] 8.809853e-01 2.380295e-01 0.1190147
[55,] 8.943004e-01 2.113992e-01 0.1056996
[56,] 8.363377e-01 3.273246e-01 0.1636623
[57,] 7.256217e-01 5.487566e-01 0.2743783
> postscript(file="/var/www/html/rcomp/tmp/1qbo71260885394.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/28urf1260885394.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/3ivcj1260885394.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/40x171260885394.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/5wepf1260885394.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 = 70
Frequency = 1
1 2 3 4 5 6
-8.355451 -124.095113 -29.251613 -154.348538 -44.886085 -95.127607
7 8 9 10 11 12
-38.339540 63.169890 35.330419 23.562056 9.049532 2.008570
13 14 15 16 17 18
44.406742 -22.516820 -32.883455 -102.524008 -10.801368 20.279295
19 20 21 22 23 24
55.663814 -13.967403 -13.639567 38.504346 104.542651 155.102970
25 26 27 28 29 30
138.479156 106.093229 -14.754097 -87.191418 -228.042412 122.588510
31 32 33 34 35 36
159.710218 116.421560 183.673906 86.448132 120.480990 188.113807
37 38 39 40 41 42
102.948803 -97.454578 275.343716 160.297558 -14.406716 -0.738291
43 44 45 46 47 48
-304.903443 130.536857 181.213896 -279.229303 45.120574 -237.156626
49 50 51 52 53 54
-102.523642 -31.538274 196.509430 -32.499812 -274.925973 -465.388127
55 56 57 58 59 60
-4.111345 -4.406412 -688.065530 -147.792292 -125.249021 61.474480
61 62 63 64 65 66
-72.419002 -123.278361 205.118241 209.706624 10.714983 73.929296
67 68 69 70
258.312245 198.507945 127.500735 15.946069
> postscript(file="/var/www/html/rcomp/tmp/6i2hr1260885394.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 = 70
Frequency = 1
lag(myerror, k = 1) myerror
0 -8.355451 NA
1 -124.095113 -8.355451
2 -29.251613 -124.095113
3 -154.348538 -29.251613
4 -44.886085 -154.348538
5 -95.127607 -44.886085
6 -38.339540 -95.127607
7 63.169890 -38.339540
8 35.330419 63.169890
9 23.562056 35.330419
10 9.049532 23.562056
11 2.008570 9.049532
12 44.406742 2.008570
13 -22.516820 44.406742
14 -32.883455 -22.516820
15 -102.524008 -32.883455
16 -10.801368 -102.524008
17 20.279295 -10.801368
18 55.663814 20.279295
19 -13.967403 55.663814
20 -13.639567 -13.967403
21 38.504346 -13.639567
22 104.542651 38.504346
23 155.102970 104.542651
24 138.479156 155.102970
25 106.093229 138.479156
26 -14.754097 106.093229
27 -87.191418 -14.754097
28 -228.042412 -87.191418
29 122.588510 -228.042412
30 159.710218 122.588510
31 116.421560 159.710218
32 183.673906 116.421560
33 86.448132 183.673906
34 120.480990 86.448132
35 188.113807 120.480990
36 102.948803 188.113807
37 -97.454578 102.948803
38 275.343716 -97.454578
39 160.297558 275.343716
40 -14.406716 160.297558
41 -0.738291 -14.406716
42 -304.903443 -0.738291
43 130.536857 -304.903443
44 181.213896 130.536857
45 -279.229303 181.213896
46 45.120574 -279.229303
47 -237.156626 45.120574
48 -102.523642 -237.156626
49 -31.538274 -102.523642
50 196.509430 -31.538274
51 -32.499812 196.509430
52 -274.925973 -32.499812
53 -465.388127 -274.925973
54 -4.111345 -465.388127
55 -4.406412 -4.111345
56 -688.065530 -4.406412
57 -147.792292 -688.065530
58 -125.249021 -147.792292
59 61.474480 -125.249021
60 -72.419002 61.474480
61 -123.278361 -72.419002
62 205.118241 -123.278361
63 209.706624 205.118241
64 10.714983 209.706624
65 73.929296 10.714983
66 258.312245 73.929296
67 198.507945 258.312245
68 127.500735 198.507945
69 15.946069 127.500735
70 NA 15.946069
> dum1 <- dum[2:length(myerror),]
> dum1
lag(myerror, k = 1) myerror
[1,] -124.095113 -8.355451
[2,] -29.251613 -124.095113
[3,] -154.348538 -29.251613
[4,] -44.886085 -154.348538
[5,] -95.127607 -44.886085
[6,] -38.339540 -95.127607
[7,] 63.169890 -38.339540
[8,] 35.330419 63.169890
[9,] 23.562056 35.330419
[10,] 9.049532 23.562056
[11,] 2.008570 9.049532
[12,] 44.406742 2.008570
[13,] -22.516820 44.406742
[14,] -32.883455 -22.516820
[15,] -102.524008 -32.883455
[16,] -10.801368 -102.524008
[17,] 20.279295 -10.801368
[18,] 55.663814 20.279295
[19,] -13.967403 55.663814
[20,] -13.639567 -13.967403
[21,] 38.504346 -13.639567
[22,] 104.542651 38.504346
[23,] 155.102970 104.542651
[24,] 138.479156 155.102970
[25,] 106.093229 138.479156
[26,] -14.754097 106.093229
[27,] -87.191418 -14.754097
[28,] -228.042412 -87.191418
[29,] 122.588510 -228.042412
[30,] 159.710218 122.588510
[31,] 116.421560 159.710218
[32,] 183.673906 116.421560
[33,] 86.448132 183.673906
[34,] 120.480990 86.448132
[35,] 188.113807 120.480990
[36,] 102.948803 188.113807
[37,] -97.454578 102.948803
[38,] 275.343716 -97.454578
[39,] 160.297558 275.343716
[40,] -14.406716 160.297558
[41,] -0.738291 -14.406716
[42,] -304.903443 -0.738291
[43,] 130.536857 -304.903443
[44,] 181.213896 130.536857
[45,] -279.229303 181.213896
[46,] 45.120574 -279.229303
[47,] -237.156626 45.120574
[48,] -102.523642 -237.156626
[49,] -31.538274 -102.523642
[50,] 196.509430 -31.538274
[51,] -32.499812 196.509430
[52,] -274.925973 -32.499812
[53,] -465.388127 -274.925973
[54,] -4.111345 -465.388127
[55,] -4.406412 -4.111345
[56,] -688.065530 -4.406412
[57,] -147.792292 -688.065530
[58,] -125.249021 -147.792292
[59,] 61.474480 -125.249021
[60,] -72.419002 61.474480
[61,] -123.278361 -72.419002
[62,] 205.118241 -123.278361
[63,] 209.706624 205.118241
[64,] 10.714983 209.706624
[65,] 73.929296 10.714983
[66,] 258.312245 73.929296
[67,] 198.507945 258.312245
[68,] 127.500735 198.507945
[69,] 15.946069 127.500735
> z <- as.data.frame(dum1)
> z
lag(myerror, k = 1) myerror
1 -124.095113 -8.355451
2 -29.251613 -124.095113
3 -154.348538 -29.251613
4 -44.886085 -154.348538
5 -95.127607 -44.886085
6 -38.339540 -95.127607
7 63.169890 -38.339540
8 35.330419 63.169890
9 23.562056 35.330419
10 9.049532 23.562056
11 2.008570 9.049532
12 44.406742 2.008570
13 -22.516820 44.406742
14 -32.883455 -22.516820
15 -102.524008 -32.883455
16 -10.801368 -102.524008
17 20.279295 -10.801368
18 55.663814 20.279295
19 -13.967403 55.663814
20 -13.639567 -13.967403
21 38.504346 -13.639567
22 104.542651 38.504346
23 155.102970 104.542651
24 138.479156 155.102970
25 106.093229 138.479156
26 -14.754097 106.093229
27 -87.191418 -14.754097
28 -228.042412 -87.191418
29 122.588510 -228.042412
30 159.710218 122.588510
31 116.421560 159.710218
32 183.673906 116.421560
33 86.448132 183.673906
34 120.480990 86.448132
35 188.113807 120.480990
36 102.948803 188.113807
37 -97.454578 102.948803
38 275.343716 -97.454578
39 160.297558 275.343716
40 -14.406716 160.297558
41 -0.738291 -14.406716
42 -304.903443 -0.738291
43 130.536857 -304.903443
44 181.213896 130.536857
45 -279.229303 181.213896
46 45.120574 -279.229303
47 -237.156626 45.120574
48 -102.523642 -237.156626
49 -31.538274 -102.523642
50 196.509430 -31.538274
51 -32.499812 196.509430
52 -274.925973 -32.499812
53 -465.388127 -274.925973
54 -4.111345 -465.388127
55 -4.406412 -4.111345
56 -688.065530 -4.406412
57 -147.792292 -688.065530
58 -125.249021 -147.792292
59 61.474480 -125.249021
60 -72.419002 61.474480
61 -123.278361 -72.419002
62 205.118241 -123.278361
63 209.706624 205.118241
64 10.714983 209.706624
65 73.929296 10.714983
66 258.312245 73.929296
67 198.507945 258.312245
68 127.500735 198.507945
69 15.946069 127.500735
> 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/727l81260885394.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/88cek1260885394.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/99i031260885394.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/1070cj1260885394.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/11166y1260885394.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/12p5rh1260885394.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/13rwlr1260885394.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/14jt3e1260885394.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/15fqfr1260885394.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/16t3qh1260885394.tab")
+ }
>
> try(system("convert tmp/1qbo71260885394.ps tmp/1qbo71260885394.png",intern=TRUE))
character(0)
> try(system("convert tmp/28urf1260885394.ps tmp/28urf1260885394.png",intern=TRUE))
character(0)
> try(system("convert tmp/3ivcj1260885394.ps tmp/3ivcj1260885394.png",intern=TRUE))
character(0)
> try(system("convert tmp/40x171260885394.ps tmp/40x171260885394.png",intern=TRUE))
character(0)
> try(system("convert tmp/5wepf1260885394.ps tmp/5wepf1260885394.png",intern=TRUE))
character(0)
> try(system("convert tmp/6i2hr1260885394.ps tmp/6i2hr1260885394.png",intern=TRUE))
character(0)
> try(system("convert tmp/727l81260885394.ps tmp/727l81260885394.png",intern=TRUE))
character(0)
> try(system("convert tmp/88cek1260885394.ps tmp/88cek1260885394.png",intern=TRUE))
character(0)
> try(system("convert tmp/99i031260885394.ps tmp/99i031260885394.png",intern=TRUE))
character(0)
> try(system("convert tmp/1070cj1260885394.ps tmp/1070cj1260885394.png",intern=TRUE))
character(0)
>
>
> proc.time()
user system elapsed
2.532 1.539 3.641