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(269285,8.2,269829,8,270911,7.5,266844,6.8,271244,6.5,269907,6.6,271296,7.6,270157,8,271322,8.1,267179,7.7,264101,7.5,265518,7.6,269419,7.8,268714,7.8,272482,7.8,268351,7.5,268175,7.5,270674,7.1,272764,7.5,272599,7.5,270333,7.6,270846,7.7,270491,7.7,269160,7.9,274027,8.1,273784,8.2,276663,8.2,274525,8.2,271344,7.9,271115,7.3,270798,6.9,273911,6.6,273985,6.7,271917,6.9,273338,7,270601,7.1,273547,7.2,275363,7.1,281229,6.9,277793,7,279913,6.8,282500,6.4,280041,6.7,282166,6.6,290304,6.4,283519,6.3,287816,6.2,285226,6.5,287595,6.8,289741,6.8,289148,6.4,288301,6.1,290155,5.8,289648,6.1,288225,7.2,289351,7.3,294735,6.9,305333,6.1,309030,5.8,310215,6.2,321935,7.1,325734,7.7,320846,7.9,323023,7.7,319753,7.4,321753,7.5,320757,8,324479,8.1),dim=c(2,68),dimnames=list(c('Y','X'),1:68))
> y <- array(NA,dim=c(2,68),dimnames=list(c('Y','X'),1:68))
> 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 269285 8.2 1 0 0 0 0 0 0 0 0 0 0 1
2 269829 8.0 0 1 0 0 0 0 0 0 0 0 0 2
3 270911 7.5 0 0 1 0 0 0 0 0 0 0 0 3
4 266844 6.8 0 0 0 1 0 0 0 0 0 0 0 4
5 271244 6.5 0 0 0 0 1 0 0 0 0 0 0 5
6 269907 6.6 0 0 0 0 0 1 0 0 0 0 0 6
7 271296 7.6 0 0 0 0 0 0 1 0 0 0 0 7
8 270157 8.0 0 0 0 0 0 0 0 1 0 0 0 8
9 271322 8.1 0 0 0 0 0 0 0 0 1 0 0 9
10 267179 7.7 0 0 0 0 0 0 0 0 0 1 0 10
11 264101 7.5 0 0 0 0 0 0 0 0 0 0 1 11
12 265518 7.6 0 0 0 0 0 0 0 0 0 0 0 12
13 269419 7.8 1 0 0 0 0 0 0 0 0 0 0 13
14 268714 7.8 0 1 0 0 0 0 0 0 0 0 0 14
15 272482 7.8 0 0 1 0 0 0 0 0 0 0 0 15
16 268351 7.5 0 0 0 1 0 0 0 0 0 0 0 16
17 268175 7.5 0 0 0 0 1 0 0 0 0 0 0 17
18 270674 7.1 0 0 0 0 0 1 0 0 0 0 0 18
19 272764 7.5 0 0 0 0 0 0 1 0 0 0 0 19
20 272599 7.5 0 0 0 0 0 0 0 1 0 0 0 20
21 270333 7.6 0 0 0 0 0 0 0 0 1 0 0 21
22 270846 7.7 0 0 0 0 0 0 0 0 0 1 0 22
23 270491 7.7 0 0 0 0 0 0 0 0 0 0 1 23
24 269160 7.9 0 0 0 0 0 0 0 0 0 0 0 24
25 274027 8.1 1 0 0 0 0 0 0 0 0 0 0 25
26 273784 8.2 0 1 0 0 0 0 0 0 0 0 0 26
27 276663 8.2 0 0 1 0 0 0 0 0 0 0 0 27
28 274525 8.2 0 0 0 1 0 0 0 0 0 0 0 28
29 271344 7.9 0 0 0 0 1 0 0 0 0 0 0 29
30 271115 7.3 0 0 0 0 0 1 0 0 0 0 0 30
31 270798 6.9 0 0 0 0 0 0 1 0 0 0 0 31
32 273911 6.6 0 0 0 0 0 0 0 1 0 0 0 32
33 273985 6.7 0 0 0 0 0 0 0 0 1 0 0 33
34 271917 6.9 0 0 0 0 0 0 0 0 0 1 0 34
35 273338 7.0 0 0 0 0 0 0 0 0 0 0 1 35
36 270601 7.1 0 0 0 0 0 0 0 0 0 0 0 36
37 273547 7.2 1 0 0 0 0 0 0 0 0 0 0 37
38 275363 7.1 0 1 0 0 0 0 0 0 0 0 0 38
39 281229 6.9 0 0 1 0 0 0 0 0 0 0 0 39
40 277793 7.0 0 0 0 1 0 0 0 0 0 0 0 40
41 279913 6.8 0 0 0 0 1 0 0 0 0 0 0 41
42 282500 6.4 0 0 0 0 0 1 0 0 0 0 0 42
43 280041 6.7 0 0 0 0 0 0 1 0 0 0 0 43
44 282166 6.6 0 0 0 0 0 0 0 1 0 0 0 44
45 290304 6.4 0 0 0 0 0 0 0 0 1 0 0 45
46 283519 6.3 0 0 0 0 0 0 0 0 0 1 0 46
47 287816 6.2 0 0 0 0 0 0 0 0 0 0 1 47
48 285226 6.5 0 0 0 0 0 0 0 0 0 0 0 48
49 287595 6.8 1 0 0 0 0 0 0 0 0 0 0 49
50 289741 6.8 0 1 0 0 0 0 0 0 0 0 0 50
51 289148 6.4 0 0 1 0 0 0 0 0 0 0 0 51
52 288301 6.1 0 0 0 1 0 0 0 0 0 0 0 52
53 290155 5.8 0 0 0 0 1 0 0 0 0 0 0 53
54 289648 6.1 0 0 0 0 0 1 0 0 0 0 0 54
55 288225 7.2 0 0 0 0 0 0 1 0 0 0 0 55
56 289351 7.3 0 0 0 0 0 0 0 1 0 0 0 56
57 294735 6.9 0 0 0 0 0 0 0 0 1 0 0 57
58 305333 6.1 0 0 0 0 0 0 0 0 0 1 0 58
59 309030 5.8 0 0 0 0 0 0 0 0 0 0 1 59
60 310215 6.2 0 0 0 0 0 0 0 0 0 0 0 60
61 321935 7.1 1 0 0 0 0 0 0 0 0 0 0 61
62 325734 7.7 0 1 0 0 0 0 0 0 0 0 0 62
63 320846 7.9 0 0 1 0 0 0 0 0 0 0 0 63
64 323023 7.7 0 0 0 1 0 0 0 0 0 0 0 64
65 319753 7.4 0 0 0 0 1 0 0 0 0 0 0 65
66 321753 7.5 0 0 0 0 0 1 0 0 0 0 0 66
67 320757 8.0 0 0 0 0 0 0 1 0 0 0 0 67
68 324479 8.1 0 0 0 0 0 0 0 1 0 0 0 68
> k <- length(x[1,])
> df <- as.data.frame(x)
> (mylm <- lm(df))
Call:
lm(formula = df)
Coefficients:
(Intercept) X M1 M2 M3 M4
209353.9 5735.5 3983.8 3986.0 5357.1 3780.1
M5 M6 M7 M8 M9 M10
4568.0 5422.2 1522.5 1953.3 2057.7 1986.2
M11 t
2914.6 841.6
> (mysum <- summary(mylm))
Call:
lm(formula = df)
Residuals:
Min 1Q Median 3Q Max
-12295 -5777 -1993 7658 16538
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 209353.87 15462.82 13.539 < 2e-16 ***
X 5735.55 1969.65 2.912 0.00521 **
M1 3983.78 5502.99 0.724 0.47223
M2 3985.99 5525.54 0.721 0.47379
M3 5357.07 5482.79 0.977 0.33289
M4 3780.11 5443.44 0.694 0.49039
M5 4567.98 5439.25 0.840 0.40471
M6 5422.22 5454.43 0.994 0.32461
M7 1522.45 5462.22 0.279 0.78152
M8 1953.35 5472.59 0.357 0.72253
M9 2057.72 5680.92 0.362 0.71860
M10 1986.24 5686.14 0.349 0.72821
M11 2914.61 5696.57 0.512 0.61099
t 841.59 61.01 13.794 < 2e-16 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 8977 on 54 degrees of freedom
Multiple R-squared: 0.789, Adjusted R-squared: 0.7382
F-statistic: 15.53 on 13 and 54 DF, p-value: 7.852e-14
> 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.361270e-02 2.722541e-02 0.98638730
[2,] 3.148427e-03 6.296854e-03 0.99685157
[3,] 9.265535e-04 1.853107e-03 0.99907345
[4,] 3.231592e-04 6.463184e-04 0.99967684
[5,] 1.344999e-04 2.689997e-04 0.99986550
[6,] 1.079409e-04 2.158819e-04 0.99989206
[7,] 3.887416e-04 7.774833e-04 0.99961126
[8,] 1.682688e-04 3.365376e-04 0.99983173
[9,] 7.983503e-05 1.596701e-04 0.99992016
[10,] 3.019251e-05 6.038502e-05 0.99996981
[11,] 1.291346e-05 2.582692e-05 0.99998709
[12,] 8.662128e-06 1.732426e-05 0.99999134
[13,] 3.653584e-06 7.307168e-06 0.99999635
[14,] 1.743938e-06 3.487877e-06 0.99999826
[15,] 2.273248e-06 4.546497e-06 0.99999773
[16,] 4.137273e-06 8.274545e-06 0.99999586
[17,] 2.776907e-06 5.553815e-06 0.99999722
[18,] 8.814425e-07 1.762885e-06 0.99999912
[19,] 6.875016e-07 1.375003e-06 0.99999931
[20,] 2.183556e-07 4.367113e-07 0.99999978
[21,] 8.265754e-08 1.653151e-07 0.99999992
[22,] 2.782297e-08 5.564595e-08 0.99999997
[23,] 3.980399e-08 7.960799e-08 0.99999996
[24,] 2.890594e-08 5.781188e-08 0.99999997
[25,] 4.652853e-08 9.305706e-08 0.99999995
[26,] 4.170725e-07 8.341451e-07 0.99999958
[27,] 1.135429e-06 2.270858e-06 0.99999886
[28,] 3.412792e-05 6.825584e-05 0.99996587
[29,] 5.026614e-01 9.946771e-01 0.49733855
[30,] 5.640458e-01 8.719083e-01 0.43595416
[31,] 8.112008e-01 3.775985e-01 0.18879923
[32,] 9.856066e-01 2.878673e-02 0.01439337
[33,] 9.697948e-01 6.041034e-02 0.03020517
[34,] 9.536151e-01 9.276970e-02 0.04638485
[35,] 8.780209e-01 2.439582e-01 0.12197908
> postscript(file="/var/www/html/rcomp/tmp/185v61258805402.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/2c9d21258805402.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/3vf3h1258805402.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/4jhpa1258805402.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/5jgw41258805402.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 = 68
Frequency = 1
1 2 3 4 5 6
8074.2756 8921.5882 10658.6964 11341.9517 15833.1549 12226.7681
7 8 9 10 11 12
10938.4020 6232.7013 5878.1870 3259.2964 -441.5489 2474.9168
13 14 15 16 17 18
403.4383 -1145.3586 409.9760 -1264.9875 -3070.4485 26.9383
19 20 21 22 23 24
2880.9006 1443.4188 -2342.0956 -3172.7597 -5297.7145 -5702.8035
25 26 27 28 29 30
-6808.2820 -8468.6336 -7802.2990 -9204.9267 -12294.7236 -10778.2273
31 32 33 34 35 36
-5742.8272 -2181.6448 -3627.1592 -7612.3781 -8534.8876 -9772.4219
37 38 39 40 41 42
-12225.3457 -10679.5878 -5879.1438 -9153.3262 -7515.6777 -4330.2909
43 44 45 46 47 48
-5451.7739 -4025.7010 4313.4488 -2668.1059 432.4941 -1805.1497
49 50 51 52 53 54
-5982.1829 -4679.9798 -5191.4263 -3582.3898 -1637.1867 -5560.6829
55 56 57 58 59 60
-10234.6037 -10954.6402 -4222.3810 10193.9474 13841.6568 14805.4583
61 62 63 64 65 66
16538.0968 16051.9716 7804.1967 11863.6785 8684.8816 8415.4948
67 68
7609.9023 9485.8659
> postscript(file="/var/www/html/rcomp/tmp/6yfcg1258805402.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 = 68
Frequency = 1
lag(myerror, k = 1) myerror
0 8074.2756 NA
1 8921.5882 8074.2756
2 10658.6964 8921.5882
3 11341.9517 10658.6964
4 15833.1549 11341.9517
5 12226.7681 15833.1549
6 10938.4020 12226.7681
7 6232.7013 10938.4020
8 5878.1870 6232.7013
9 3259.2964 5878.1870
10 -441.5489 3259.2964
11 2474.9168 -441.5489
12 403.4383 2474.9168
13 -1145.3586 403.4383
14 409.9760 -1145.3586
15 -1264.9875 409.9760
16 -3070.4485 -1264.9875
17 26.9383 -3070.4485
18 2880.9006 26.9383
19 1443.4188 2880.9006
20 -2342.0956 1443.4188
21 -3172.7597 -2342.0956
22 -5297.7145 -3172.7597
23 -5702.8035 -5297.7145
24 -6808.2820 -5702.8035
25 -8468.6336 -6808.2820
26 -7802.2990 -8468.6336
27 -9204.9267 -7802.2990
28 -12294.7236 -9204.9267
29 -10778.2273 -12294.7236
30 -5742.8272 -10778.2273
31 -2181.6448 -5742.8272
32 -3627.1592 -2181.6448
33 -7612.3781 -3627.1592
34 -8534.8876 -7612.3781
35 -9772.4219 -8534.8876
36 -12225.3457 -9772.4219
37 -10679.5878 -12225.3457
38 -5879.1438 -10679.5878
39 -9153.3262 -5879.1438
40 -7515.6777 -9153.3262
41 -4330.2909 -7515.6777
42 -5451.7739 -4330.2909
43 -4025.7010 -5451.7739
44 4313.4488 -4025.7010
45 -2668.1059 4313.4488
46 432.4941 -2668.1059
47 -1805.1497 432.4941
48 -5982.1829 -1805.1497
49 -4679.9798 -5982.1829
50 -5191.4263 -4679.9798
51 -3582.3898 -5191.4263
52 -1637.1867 -3582.3898
53 -5560.6829 -1637.1867
54 -10234.6037 -5560.6829
55 -10954.6402 -10234.6037
56 -4222.3810 -10954.6402
57 10193.9474 -4222.3810
58 13841.6568 10193.9474
59 14805.4583 13841.6568
60 16538.0968 14805.4583
61 16051.9716 16538.0968
62 7804.1967 16051.9716
63 11863.6785 7804.1967
64 8684.8816 11863.6785
65 8415.4948 8684.8816
66 7609.9023 8415.4948
67 9485.8659 7609.9023
68 NA 9485.8659
> dum1 <- dum[2:length(myerror),]
> dum1
lag(myerror, k = 1) myerror
[1,] 8921.5882 8074.2756
[2,] 10658.6964 8921.5882
[3,] 11341.9517 10658.6964
[4,] 15833.1549 11341.9517
[5,] 12226.7681 15833.1549
[6,] 10938.4020 12226.7681
[7,] 6232.7013 10938.4020
[8,] 5878.1870 6232.7013
[9,] 3259.2964 5878.1870
[10,] -441.5489 3259.2964
[11,] 2474.9168 -441.5489
[12,] 403.4383 2474.9168
[13,] -1145.3586 403.4383
[14,] 409.9760 -1145.3586
[15,] -1264.9875 409.9760
[16,] -3070.4485 -1264.9875
[17,] 26.9383 -3070.4485
[18,] 2880.9006 26.9383
[19,] 1443.4188 2880.9006
[20,] -2342.0956 1443.4188
[21,] -3172.7597 -2342.0956
[22,] -5297.7145 -3172.7597
[23,] -5702.8035 -5297.7145
[24,] -6808.2820 -5702.8035
[25,] -8468.6336 -6808.2820
[26,] -7802.2990 -8468.6336
[27,] -9204.9267 -7802.2990
[28,] -12294.7236 -9204.9267
[29,] -10778.2273 -12294.7236
[30,] -5742.8272 -10778.2273
[31,] -2181.6448 -5742.8272
[32,] -3627.1592 -2181.6448
[33,] -7612.3781 -3627.1592
[34,] -8534.8876 -7612.3781
[35,] -9772.4219 -8534.8876
[36,] -12225.3457 -9772.4219
[37,] -10679.5878 -12225.3457
[38,] -5879.1438 -10679.5878
[39,] -9153.3262 -5879.1438
[40,] -7515.6777 -9153.3262
[41,] -4330.2909 -7515.6777
[42,] -5451.7739 -4330.2909
[43,] -4025.7010 -5451.7739
[44,] 4313.4488 -4025.7010
[45,] -2668.1059 4313.4488
[46,] 432.4941 -2668.1059
[47,] -1805.1497 432.4941
[48,] -5982.1829 -1805.1497
[49,] -4679.9798 -5982.1829
[50,] -5191.4263 -4679.9798
[51,] -3582.3898 -5191.4263
[52,] -1637.1867 -3582.3898
[53,] -5560.6829 -1637.1867
[54,] -10234.6037 -5560.6829
[55,] -10954.6402 -10234.6037
[56,] -4222.3810 -10954.6402
[57,] 10193.9474 -4222.3810
[58,] 13841.6568 10193.9474
[59,] 14805.4583 13841.6568
[60,] 16538.0968 14805.4583
[61,] 16051.9716 16538.0968
[62,] 7804.1967 16051.9716
[63,] 11863.6785 7804.1967
[64,] 8684.8816 11863.6785
[65,] 8415.4948 8684.8816
[66,] 7609.9023 8415.4948
[67,] 9485.8659 7609.9023
> z <- as.data.frame(dum1)
> z
lag(myerror, k = 1) myerror
1 8921.5882 8074.2756
2 10658.6964 8921.5882
3 11341.9517 10658.6964
4 15833.1549 11341.9517
5 12226.7681 15833.1549
6 10938.4020 12226.7681
7 6232.7013 10938.4020
8 5878.1870 6232.7013
9 3259.2964 5878.1870
10 -441.5489 3259.2964
11 2474.9168 -441.5489
12 403.4383 2474.9168
13 -1145.3586 403.4383
14 409.9760 -1145.3586
15 -1264.9875 409.9760
16 -3070.4485 -1264.9875
17 26.9383 -3070.4485
18 2880.9006 26.9383
19 1443.4188 2880.9006
20 -2342.0956 1443.4188
21 -3172.7597 -2342.0956
22 -5297.7145 -3172.7597
23 -5702.8035 -5297.7145
24 -6808.2820 -5702.8035
25 -8468.6336 -6808.2820
26 -7802.2990 -8468.6336
27 -9204.9267 -7802.2990
28 -12294.7236 -9204.9267
29 -10778.2273 -12294.7236
30 -5742.8272 -10778.2273
31 -2181.6448 -5742.8272
32 -3627.1592 -2181.6448
33 -7612.3781 -3627.1592
34 -8534.8876 -7612.3781
35 -9772.4219 -8534.8876
36 -12225.3457 -9772.4219
37 -10679.5878 -12225.3457
38 -5879.1438 -10679.5878
39 -9153.3262 -5879.1438
40 -7515.6777 -9153.3262
41 -4330.2909 -7515.6777
42 -5451.7739 -4330.2909
43 -4025.7010 -5451.7739
44 4313.4488 -4025.7010
45 -2668.1059 4313.4488
46 432.4941 -2668.1059
47 -1805.1497 432.4941
48 -5982.1829 -1805.1497
49 -4679.9798 -5982.1829
50 -5191.4263 -4679.9798
51 -3582.3898 -5191.4263
52 -1637.1867 -3582.3898
53 -5560.6829 -1637.1867
54 -10234.6037 -5560.6829
55 -10954.6402 -10234.6037
56 -4222.3810 -10954.6402
57 10193.9474 -4222.3810
58 13841.6568 10193.9474
59 14805.4583 13841.6568
60 16538.0968 14805.4583
61 16051.9716 16538.0968
62 7804.1967 16051.9716
63 11863.6785 7804.1967
64 8684.8816 11863.6785
65 8415.4948 8684.8816
66 7609.9023 8415.4948
67 9485.8659 7609.9023
> 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/7no0o1258805402.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/8dooo1258805402.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/9xt671258805402.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/1009ix1258805402.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/11bo9u1258805402.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/12aj9r1258805402.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/13k51z1258805402.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/1462o41258805402.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/150h0a1258805402.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/1666781258805402.tab")
+ }
>
> system("convert tmp/185v61258805402.ps tmp/185v61258805402.png")
> system("convert tmp/2c9d21258805402.ps tmp/2c9d21258805402.png")
> system("convert tmp/3vf3h1258805402.ps tmp/3vf3h1258805402.png")
> system("convert tmp/4jhpa1258805402.ps tmp/4jhpa1258805402.png")
> system("convert tmp/5jgw41258805402.ps tmp/5jgw41258805402.png")
> system("convert tmp/6yfcg1258805402.ps tmp/6yfcg1258805402.png")
> system("convert tmp/7no0o1258805402.ps tmp/7no0o1258805402.png")
> system("convert tmp/8dooo1258805402.ps tmp/8dooo1258805402.png")
> system("convert tmp/9xt671258805402.ps tmp/9xt671258805402.png")
> system("convert tmp/1009ix1258805402.ps tmp/1009ix1258805402.png")
>
>
> proc.time()
user system elapsed
2.527 1.613 3.809