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(467,98.8,460,100.5,448,110.4,443,96.4,436,101.9,431,106.2,484,81,510,94.7,513,101,503,109.4,471,102.3,471,90.7,476,96.2,475,96.1,470,106,461,103.1,455,102,456,104.7,517,86,525,92.1,523,106.9,519,112.6,509,101.7,512,92,519,97.4,517,97,510,105.4,509,102.7,501,98.1,507,104.5,569,87.4,580,89.9,578,109.8,565,111.7,547,98.6,555,96.9,562,95.1,561,97,555,112.7,544,102.9,537,97.4,543,111.4,594,87.4,611,96.8,613,114.1,611,110.3,594,103.9,595,101.6,591,94.6,589,95.9,584,104.7,573,102.8,567,98.1,569,113.9,621,80.9,629,95.7,628,113.2,612,105.9,595,108.8,597,102.3,593,99,590,100.7,580,115.5,574,100.7,573,109.9,573,114.6,620,85.4,626,100.5,620,114.8,588,116.5,566,112.9,557,102),dim=c(2,72),dimnames=list(c('Y','X'),1:72))
> y <- array(NA,dim=c(2,72),dimnames=list(c('Y','X'),1:72))
> 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 467 98.8 1 0 0 0 0 0 0 0 0 0 0 1
2 460 100.5 0 1 0 0 0 0 0 0 0 0 0 2
3 448 110.4 0 0 1 0 0 0 0 0 0 0 0 3
4 443 96.4 0 0 0 1 0 0 0 0 0 0 0 4
5 436 101.9 0 0 0 0 1 0 0 0 0 0 0 5
6 431 106.2 0 0 0 0 0 1 0 0 0 0 0 6
7 484 81.0 0 0 0 0 0 0 1 0 0 0 0 7
8 510 94.7 0 0 0 0 0 0 0 1 0 0 0 8
9 513 101.0 0 0 0 0 0 0 0 0 1 0 0 9
10 503 109.4 0 0 0 0 0 0 0 0 0 1 0 10
11 471 102.3 0 0 0 0 0 0 0 0 0 0 1 11
12 471 90.7 0 0 0 0 0 0 0 0 0 0 0 12
13 476 96.2 1 0 0 0 0 0 0 0 0 0 0 13
14 475 96.1 0 1 0 0 0 0 0 0 0 0 0 14
15 470 106.0 0 0 1 0 0 0 0 0 0 0 0 15
16 461 103.1 0 0 0 1 0 0 0 0 0 0 0 16
17 455 102.0 0 0 0 0 1 0 0 0 0 0 0 17
18 456 104.7 0 0 0 0 0 1 0 0 0 0 0 18
19 517 86.0 0 0 0 0 0 0 1 0 0 0 0 19
20 525 92.1 0 0 0 0 0 0 0 1 0 0 0 20
21 523 106.9 0 0 0 0 0 0 0 0 1 0 0 21
22 519 112.6 0 0 0 0 0 0 0 0 0 1 0 22
23 509 101.7 0 0 0 0 0 0 0 0 0 0 1 23
24 512 92.0 0 0 0 0 0 0 0 0 0 0 0 24
25 519 97.4 1 0 0 0 0 0 0 0 0 0 0 25
26 517 97.0 0 1 0 0 0 0 0 0 0 0 0 26
27 510 105.4 0 0 1 0 0 0 0 0 0 0 0 27
28 509 102.7 0 0 0 1 0 0 0 0 0 0 0 28
29 501 98.1 0 0 0 0 1 0 0 0 0 0 0 29
30 507 104.5 0 0 0 0 0 1 0 0 0 0 0 30
31 569 87.4 0 0 0 0 0 0 1 0 0 0 0 31
32 580 89.9 0 0 0 0 0 0 0 1 0 0 0 32
33 578 109.8 0 0 0 0 0 0 0 0 1 0 0 33
34 565 111.7 0 0 0 0 0 0 0 0 0 1 0 34
35 547 98.6 0 0 0 0 0 0 0 0 0 0 1 35
36 555 96.9 0 0 0 0 0 0 0 0 0 0 0 36
37 562 95.1 1 0 0 0 0 0 0 0 0 0 0 37
38 561 97.0 0 1 0 0 0 0 0 0 0 0 0 38
39 555 112.7 0 0 1 0 0 0 0 0 0 0 0 39
40 544 102.9 0 0 0 1 0 0 0 0 0 0 0 40
41 537 97.4 0 0 0 0 1 0 0 0 0 0 0 41
42 543 111.4 0 0 0 0 0 1 0 0 0 0 0 42
43 594 87.4 0 0 0 0 0 0 1 0 0 0 0 43
44 611 96.8 0 0 0 0 0 0 0 1 0 0 0 44
45 613 114.1 0 0 0 0 0 0 0 0 1 0 0 45
46 611 110.3 0 0 0 0 0 0 0 0 0 1 0 46
47 594 103.9 0 0 0 0 0 0 0 0 0 0 1 47
48 595 101.6 0 0 0 0 0 0 0 0 0 0 0 48
49 591 94.6 1 0 0 0 0 0 0 0 0 0 0 49
50 589 95.9 0 1 0 0 0 0 0 0 0 0 0 50
51 584 104.7 0 0 1 0 0 0 0 0 0 0 0 51
52 573 102.8 0 0 0 1 0 0 0 0 0 0 0 52
53 567 98.1 0 0 0 0 1 0 0 0 0 0 0 53
54 569 113.9 0 0 0 0 0 1 0 0 0 0 0 54
55 621 80.9 0 0 0 0 0 0 1 0 0 0 0 55
56 629 95.7 0 0 0 0 0 0 0 1 0 0 0 56
57 628 113.2 0 0 0 0 0 0 0 0 1 0 0 57
58 612 105.9 0 0 0 0 0 0 0 0 0 1 0 58
59 595 108.8 0 0 0 0 0 0 0 0 0 0 1 59
60 597 102.3 0 0 0 0 0 0 0 0 0 0 0 60
61 593 99.0 1 0 0 0 0 0 0 0 0 0 0 61
62 590 100.7 0 1 0 0 0 0 0 0 0 0 0 62
63 580 115.5 0 0 1 0 0 0 0 0 0 0 0 63
64 574 100.7 0 0 0 1 0 0 0 0 0 0 0 64
65 573 109.9 0 0 0 0 1 0 0 0 0 0 0 65
66 573 114.6 0 0 0 0 0 1 0 0 0 0 0 66
67 620 85.4 0 0 0 0 0 0 1 0 0 0 0 67
68 626 100.5 0 0 0 0 0 0 0 1 0 0 0 68
69 620 114.8 0 0 0 0 0 0 0 0 1 0 0 69
70 588 116.5 0 0 0 0 0 0 0 0 0 1 0 70
71 566 112.9 0 0 0 0 0 0 0 0 0 0 1 71
72 557 102.0 0 0 0 0 0 0 0 0 0 0 0 72
> k <- length(x[1,])
> df <- as.data.frame(x)
> (mylm <- lm(df))
Call:
lm(formula = df)
Coefficients:
(Intercept) X M1 M2 M3 M4
586.780 -1.427 12.049 8.446 14.609 -5.907
M5 M6 M7 M8 M9 M10
-14.413 -3.744 13.200 38.126 56.163 42.512
M11 t
11.707 2.387
> (mysum <- summary(mylm))
Call:
lm(formula = df)
Residuals:
Min 1Q Median 3Q Max
-56.157 -10.413 0.898 10.961 38.572
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 586.7804 66.4216 8.834 2.49e-12 ***
X -1.4267 0.7034 -2.028 0.047119 *
M1 12.0495 10.8688 1.109 0.272164
M2 8.4458 10.8900 0.776 0.441163
M3 14.6086 13.9124 1.050 0.298054
M4 -5.9073 11.3111 -0.522 0.603480
M5 -14.4135 11.2456 -1.282 0.205048
M6 -3.7445 13.8123 -0.271 0.787278
M7 13.1998 13.9046 0.949 0.346402
M8 38.1263 10.9260 3.490 0.000931 ***
M9 56.1631 14.0066 4.010 0.000176 ***
M10 42.5116 14.4645 2.939 0.004720 **
M11 11.7075 11.9341 0.981 0.330662
t 2.3875 0.1267 18.844 < 2e-16 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 18.71 on 58 degrees of freedom
Multiple R-squared: 0.9063, Adjusted R-squared: 0.8853
F-statistic: 43.17 on 13 and 58 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,] 0.0109086222 0.021817244 0.9890914
[2,] 0.0052491125 0.010498225 0.9947509
[3,] 0.0046317893 0.009263579 0.9953682
[4,] 0.0014575409 0.002915082 0.9985425
[5,] 0.0016965154 0.003393031 0.9983035
[6,] 0.0005574717 0.001114943 0.9994425
[7,] 0.0018009239 0.003601848 0.9981991
[8,] 0.0039234479 0.007846896 0.9960766
[9,] 0.0041039264 0.008207853 0.9958961
[10,] 0.0042861138 0.008572228 0.9957139
[11,] 0.0048265985 0.009653197 0.9951734
[12,] 0.0067080988 0.013416198 0.9932919
[13,] 0.0079681579 0.015936316 0.9920318
[14,] 0.0183020346 0.036604069 0.9816980
[15,] 0.0304236075 0.060847215 0.9695764
[16,] 0.0352910818 0.070582164 0.9647089
[17,] 0.0349272256 0.069854451 0.9650728
[18,] 0.0257304227 0.051460845 0.9742696
[19,] 0.0313800724 0.062760145 0.9686199
[20,] 0.0296818421 0.059363684 0.9703182
[21,] 0.0265457225 0.053091445 0.9734543
[22,] 0.0244250439 0.048850088 0.9755750
[23,] 0.0194568363 0.038913673 0.9805432
[24,] 0.0181115950 0.036223190 0.9818884
[25,] 0.0413472168 0.082694434 0.9586528
[26,] 0.0761551610 0.152310322 0.9238448
[27,] 0.1177621110 0.235524222 0.8822379
[28,] 0.1817382493 0.363476499 0.8182618
[29,] 0.2842737985 0.568547597 0.7157262
[30,] 0.2404157618 0.480831524 0.7595842
[31,] 0.1962038910 0.392407782 0.8037961
[32,] 0.1421836452 0.284367290 0.8578164
[33,] 0.1120677382 0.224135476 0.8879323
[34,] 0.0845419156 0.169083831 0.9154581
[35,] 0.0510994930 0.102198986 0.9489005
[36,] 0.0501014662 0.100202932 0.9498985
[37,] 0.0462801262 0.092560252 0.9537199
[38,] 0.0728243324 0.145648665 0.9271757
[39,] 0.0896965841 0.179393168 0.9103034
> postscript(file="/var/www/html/rcomp/tmp/1jpax1260636589.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/227cv1260636589.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/3vf3w1260636589.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/4htln1260636589.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/5omuc1260636589.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 = 72
Frequency = 1
1 2 3 4 5 6
6.7401392 3.3817146 -3.0443245 -9.8896187 -2.9241202 -14.8457825
7 8 9 10 11 12
-17.1302463 1.1014088 -7.3346216 5.9135900 -7.7993202 -15.0290057
13 14 15 16 17 18
-16.6191692 -16.5456460 -15.9716852 -10.9806576 -12.4313504 -20.6357257
19 20 21 22 23 24
-5.6466678 -16.2578996 -17.5670170 -2.1708837 0.6947627 -0.8242011
25 26 27 28 29 30
-0.5570342 -1.9115197 -5.4776023 7.7987644 -0.6453632 1.4290354
31 32 33 34 35 36
19.7008064 6.9534702 12.9205006 13.8951905 5.6221064 20.5167078
37 38 39 40 41 42
10.5116661 13.4385806 21.2873762 14.4342037 5.7060501 18.6233356
43 44 45 46 47 48
16.0509066 19.1477704 25.4053921 29.2479168 31.5336936 38.5722776
49 50 51 52 53 54
10.1484185 11.2193156 10.2239113 14.6416344 8.0548373 19.5401750
55 56 57 58 59 60
5.1274852 6.9285054 10.4714663 -4.6794438 10.8746025 12.9210649
61 62 63 64 65 66
-10.2240204 -9.5824450 -7.0176755 -16.0043262 2.2399463 -4.1110378
67 68 69 70 71 72
-18.1022842 -17.8732552 -23.8957204 -42.2063697 -40.9258451 -56.1568436
> postscript(file="/var/www/html/rcomp/tmp/6v12o1260636589.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 = 72
Frequency = 1
lag(myerror, k = 1) myerror
0 6.7401392 NA
1 3.3817146 6.7401392
2 -3.0443245 3.3817146
3 -9.8896187 -3.0443245
4 -2.9241202 -9.8896187
5 -14.8457825 -2.9241202
6 -17.1302463 -14.8457825
7 1.1014088 -17.1302463
8 -7.3346216 1.1014088
9 5.9135900 -7.3346216
10 -7.7993202 5.9135900
11 -15.0290057 -7.7993202
12 -16.6191692 -15.0290057
13 -16.5456460 -16.6191692
14 -15.9716852 -16.5456460
15 -10.9806576 -15.9716852
16 -12.4313504 -10.9806576
17 -20.6357257 -12.4313504
18 -5.6466678 -20.6357257
19 -16.2578996 -5.6466678
20 -17.5670170 -16.2578996
21 -2.1708837 -17.5670170
22 0.6947627 -2.1708837
23 -0.8242011 0.6947627
24 -0.5570342 -0.8242011
25 -1.9115197 -0.5570342
26 -5.4776023 -1.9115197
27 7.7987644 -5.4776023
28 -0.6453632 7.7987644
29 1.4290354 -0.6453632
30 19.7008064 1.4290354
31 6.9534702 19.7008064
32 12.9205006 6.9534702
33 13.8951905 12.9205006
34 5.6221064 13.8951905
35 20.5167078 5.6221064
36 10.5116661 20.5167078
37 13.4385806 10.5116661
38 21.2873762 13.4385806
39 14.4342037 21.2873762
40 5.7060501 14.4342037
41 18.6233356 5.7060501
42 16.0509066 18.6233356
43 19.1477704 16.0509066
44 25.4053921 19.1477704
45 29.2479168 25.4053921
46 31.5336936 29.2479168
47 38.5722776 31.5336936
48 10.1484185 38.5722776
49 11.2193156 10.1484185
50 10.2239113 11.2193156
51 14.6416344 10.2239113
52 8.0548373 14.6416344
53 19.5401750 8.0548373
54 5.1274852 19.5401750
55 6.9285054 5.1274852
56 10.4714663 6.9285054
57 -4.6794438 10.4714663
58 10.8746025 -4.6794438
59 12.9210649 10.8746025
60 -10.2240204 12.9210649
61 -9.5824450 -10.2240204
62 -7.0176755 -9.5824450
63 -16.0043262 -7.0176755
64 2.2399463 -16.0043262
65 -4.1110378 2.2399463
66 -18.1022842 -4.1110378
67 -17.8732552 -18.1022842
68 -23.8957204 -17.8732552
69 -42.2063697 -23.8957204
70 -40.9258451 -42.2063697
71 -56.1568436 -40.9258451
72 NA -56.1568436
> dum1 <- dum[2:length(myerror),]
> dum1
lag(myerror, k = 1) myerror
[1,] 3.3817146 6.7401392
[2,] -3.0443245 3.3817146
[3,] -9.8896187 -3.0443245
[4,] -2.9241202 -9.8896187
[5,] -14.8457825 -2.9241202
[6,] -17.1302463 -14.8457825
[7,] 1.1014088 -17.1302463
[8,] -7.3346216 1.1014088
[9,] 5.9135900 -7.3346216
[10,] -7.7993202 5.9135900
[11,] -15.0290057 -7.7993202
[12,] -16.6191692 -15.0290057
[13,] -16.5456460 -16.6191692
[14,] -15.9716852 -16.5456460
[15,] -10.9806576 -15.9716852
[16,] -12.4313504 -10.9806576
[17,] -20.6357257 -12.4313504
[18,] -5.6466678 -20.6357257
[19,] -16.2578996 -5.6466678
[20,] -17.5670170 -16.2578996
[21,] -2.1708837 -17.5670170
[22,] 0.6947627 -2.1708837
[23,] -0.8242011 0.6947627
[24,] -0.5570342 -0.8242011
[25,] -1.9115197 -0.5570342
[26,] -5.4776023 -1.9115197
[27,] 7.7987644 -5.4776023
[28,] -0.6453632 7.7987644
[29,] 1.4290354 -0.6453632
[30,] 19.7008064 1.4290354
[31,] 6.9534702 19.7008064
[32,] 12.9205006 6.9534702
[33,] 13.8951905 12.9205006
[34,] 5.6221064 13.8951905
[35,] 20.5167078 5.6221064
[36,] 10.5116661 20.5167078
[37,] 13.4385806 10.5116661
[38,] 21.2873762 13.4385806
[39,] 14.4342037 21.2873762
[40,] 5.7060501 14.4342037
[41,] 18.6233356 5.7060501
[42,] 16.0509066 18.6233356
[43,] 19.1477704 16.0509066
[44,] 25.4053921 19.1477704
[45,] 29.2479168 25.4053921
[46,] 31.5336936 29.2479168
[47,] 38.5722776 31.5336936
[48,] 10.1484185 38.5722776
[49,] 11.2193156 10.1484185
[50,] 10.2239113 11.2193156
[51,] 14.6416344 10.2239113
[52,] 8.0548373 14.6416344
[53,] 19.5401750 8.0548373
[54,] 5.1274852 19.5401750
[55,] 6.9285054 5.1274852
[56,] 10.4714663 6.9285054
[57,] -4.6794438 10.4714663
[58,] 10.8746025 -4.6794438
[59,] 12.9210649 10.8746025
[60,] -10.2240204 12.9210649
[61,] -9.5824450 -10.2240204
[62,] -7.0176755 -9.5824450
[63,] -16.0043262 -7.0176755
[64,] 2.2399463 -16.0043262
[65,] -4.1110378 2.2399463
[66,] -18.1022842 -4.1110378
[67,] -17.8732552 -18.1022842
[68,] -23.8957204 -17.8732552
[69,] -42.2063697 -23.8957204
[70,] -40.9258451 -42.2063697
[71,] -56.1568436 -40.9258451
> z <- as.data.frame(dum1)
> z
lag(myerror, k = 1) myerror
1 3.3817146 6.7401392
2 -3.0443245 3.3817146
3 -9.8896187 -3.0443245
4 -2.9241202 -9.8896187
5 -14.8457825 -2.9241202
6 -17.1302463 -14.8457825
7 1.1014088 -17.1302463
8 -7.3346216 1.1014088
9 5.9135900 -7.3346216
10 -7.7993202 5.9135900
11 -15.0290057 -7.7993202
12 -16.6191692 -15.0290057
13 -16.5456460 -16.6191692
14 -15.9716852 -16.5456460
15 -10.9806576 -15.9716852
16 -12.4313504 -10.9806576
17 -20.6357257 -12.4313504
18 -5.6466678 -20.6357257
19 -16.2578996 -5.6466678
20 -17.5670170 -16.2578996
21 -2.1708837 -17.5670170
22 0.6947627 -2.1708837
23 -0.8242011 0.6947627
24 -0.5570342 -0.8242011
25 -1.9115197 -0.5570342
26 -5.4776023 -1.9115197
27 7.7987644 -5.4776023
28 -0.6453632 7.7987644
29 1.4290354 -0.6453632
30 19.7008064 1.4290354
31 6.9534702 19.7008064
32 12.9205006 6.9534702
33 13.8951905 12.9205006
34 5.6221064 13.8951905
35 20.5167078 5.6221064
36 10.5116661 20.5167078
37 13.4385806 10.5116661
38 21.2873762 13.4385806
39 14.4342037 21.2873762
40 5.7060501 14.4342037
41 18.6233356 5.7060501
42 16.0509066 18.6233356
43 19.1477704 16.0509066
44 25.4053921 19.1477704
45 29.2479168 25.4053921
46 31.5336936 29.2479168
47 38.5722776 31.5336936
48 10.1484185 38.5722776
49 11.2193156 10.1484185
50 10.2239113 11.2193156
51 14.6416344 10.2239113
52 8.0548373 14.6416344
53 19.5401750 8.0548373
54 5.1274852 19.5401750
55 6.9285054 5.1274852
56 10.4714663 6.9285054
57 -4.6794438 10.4714663
58 10.8746025 -4.6794438
59 12.9210649 10.8746025
60 -10.2240204 12.9210649
61 -9.5824450 -10.2240204
62 -7.0176755 -9.5824450
63 -16.0043262 -7.0176755
64 2.2399463 -16.0043262
65 -4.1110378 2.2399463
66 -18.1022842 -4.1110378
67 -17.8732552 -18.1022842
68 -23.8957204 -17.8732552
69 -42.2063697 -23.8957204
70 -40.9258451 -42.2063697
71 -56.1568436 -40.9258451
> 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/7h4rg1260636589.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/8z14u1260636589.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/9oqnn1260636589.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/10t7o01260636589.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/112wjv1260636589.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/12057l1260636589.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/13m6161260636589.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/14drqc1260636589.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/153xdt1260636589.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/167ume1260636589.tab")
+ }
>
> try(system("convert tmp/1jpax1260636589.ps tmp/1jpax1260636589.png",intern=TRUE))
character(0)
> try(system("convert tmp/227cv1260636589.ps tmp/227cv1260636589.png",intern=TRUE))
character(0)
> try(system("convert tmp/3vf3w1260636589.ps tmp/3vf3w1260636589.png",intern=TRUE))
character(0)
> try(system("convert tmp/4htln1260636589.ps tmp/4htln1260636589.png",intern=TRUE))
character(0)
> try(system("convert tmp/5omuc1260636589.ps tmp/5omuc1260636589.png",intern=TRUE))
character(0)
> try(system("convert tmp/6v12o1260636589.ps tmp/6v12o1260636589.png",intern=TRUE))
character(0)
> try(system("convert tmp/7h4rg1260636589.ps tmp/7h4rg1260636589.png",intern=TRUE))
character(0)
> try(system("convert tmp/8z14u1260636589.ps tmp/8z14u1260636589.png",intern=TRUE))
character(0)
> try(system("convert tmp/9oqnn1260636589.ps tmp/9oqnn1260636589.png",intern=TRUE))
character(0)
> try(system("convert tmp/10t7o01260636589.ps tmp/10t7o01260636589.png",intern=TRUE))
character(0)
>
>
> proc.time()
user system elapsed
2.510 1.545 3.191