R version 2.8.0 (2008-10-20)
Copyright (C) 2008 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
+ ,101.0
+ ,0
+ ,460
+ ,98.7
+ ,0
+ ,448
+ ,105.1
+ ,0
+ ,443
+ ,98.4
+ ,0
+ ,436
+ ,101.7
+ ,0
+ ,431
+ ,102.9
+ ,0
+ ,484
+ ,92.2
+ ,0
+ ,510
+ ,94.9
+ ,0
+ ,513
+ ,92.8
+ ,0
+ ,503
+ ,98.5
+ ,0
+ ,471
+ ,94.3
+ ,0
+ ,471
+ ,87.4
+ ,0
+ ,476
+ ,103.4
+ ,0
+ ,475
+ ,101.2
+ ,0
+ ,470
+ ,109.6
+ ,0
+ ,461
+ ,111.9
+ ,0
+ ,455
+ ,108.9
+ ,0
+ ,456
+ ,105.6
+ ,0
+ ,517
+ ,107.8
+ ,0
+ ,525
+ ,97.5
+ ,0
+ ,523
+ ,102.4
+ ,0
+ ,519
+ ,105.6
+ ,0
+ ,509
+ ,99.8
+ ,0
+ ,512
+ ,96.2
+ ,0
+ ,519
+ ,113.1
+ ,0
+ ,517
+ ,107.4
+ ,0
+ ,510
+ ,116.8
+ ,0
+ ,509
+ ,112.9
+ ,0
+ ,501
+ ,105.3
+ ,0
+ ,507
+ ,109.3
+ ,0
+ ,569
+ ,107.9
+ ,0
+ ,580
+ ,101.1
+ ,0
+ ,578
+ ,114.7
+ ,0
+ ,565
+ ,116.2
+ ,0
+ ,547
+ ,108.4
+ ,0
+ ,555
+ ,113.4
+ ,0
+ ,562
+ ,108.7
+ ,0
+ ,561
+ ,112.6
+ ,0
+ ,555
+ ,124.2
+ ,1
+ ,544
+ ,114.9
+ ,1
+ ,537
+ ,110.5
+ ,1
+ ,543
+ ,121.5
+ ,1
+ ,594
+ ,118.1
+ ,1
+ ,611
+ ,111.7
+ ,1
+ ,613
+ ,132.7
+ ,1
+ ,611
+ ,119.0
+ ,1
+ ,594
+ ,116.7
+ ,1
+ ,595
+ ,120.1
+ ,1
+ ,591
+ ,113.4
+ ,1
+ ,589
+ ,106.6
+ ,1
+ ,584
+ ,116.3
+ ,1
+ ,573
+ ,112.6
+ ,1
+ ,567
+ ,111.6
+ ,1
+ ,569
+ ,125.1
+ ,1
+ ,621
+ ,110.7
+ ,1
+ ,629
+ ,109.6
+ ,1
+ ,628
+ ,114.2
+ ,1
+ ,612
+ ,113.4
+ ,1
+ ,595
+ ,116.0
+ ,1
+ ,597
+ ,109.6
+ ,1
+ ,593
+ ,117.8
+ ,1
+ ,590
+ ,115.8
+ ,1
+ ,580
+ ,125.3
+ ,1
+ ,574
+ ,113.0
+ ,1
+ ,573
+ ,120.5
+ ,1
+ ,573
+ ,116.6
+ ,1
+ ,620
+ ,111.8
+ ,1
+ ,626
+ ,115.2
+ ,1
+ ,620
+ ,118.6
+ ,1
+ ,588
+ ,122.4
+ ,1
+ ,566
+ ,116.4
+ ,1
+ ,557
+ ,114.5
+ ,1
+ ,561
+ ,119.8
+ ,1
+ ,549
+ ,115.8
+ ,1
+ ,532
+ ,127.8
+ ,1
+ ,526
+ ,118.8
+ ,1
+ ,511
+ ,119.7
+ ,1
+ ,499
+ ,118.6
+ ,1
+ ,555
+ ,120.8
+ ,1
+ ,565
+ ,115.9
+ ,1
+ ,542
+ ,109.7
+ ,1
+ ,527
+ ,114.8
+ ,1
+ ,510
+ ,116.2
+ ,1
+ ,514
+ ,112.2
+ ,1)
+ ,dim=c(3
+ ,84)
+ ,dimnames=list(c('Y'
+ ,'X'
+ ,'DUM')
+ ,1:84))
> y <- array(NA,dim=c(3,84),dimnames=list(c('Y','X','DUM'),1:84))
> 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 DUM M1 M2 M3 M4 M5 M6 M7 M8 M9 M10 M11 t
1 467 101.0 0 1 0 0 0 0 0 0 0 0 0 0 1
2 460 98.7 0 0 1 0 0 0 0 0 0 0 0 0 2
3 448 105.1 0 0 0 1 0 0 0 0 0 0 0 0 3
4 443 98.4 0 0 0 0 1 0 0 0 0 0 0 0 4
5 436 101.7 0 0 0 0 0 1 0 0 0 0 0 0 5
6 431 102.9 0 0 0 0 0 0 1 0 0 0 0 0 6
7 484 92.2 0 0 0 0 0 0 0 1 0 0 0 0 7
8 510 94.9 0 0 0 0 0 0 0 0 1 0 0 0 8
9 513 92.8 0 0 0 0 0 0 0 0 0 1 0 0 9
10 503 98.5 0 0 0 0 0 0 0 0 0 0 1 0 10
11 471 94.3 0 0 0 0 0 0 0 0 0 0 0 1 11
12 471 87.4 0 0 0 0 0 0 0 0 0 0 0 0 12
13 476 103.4 0 1 0 0 0 0 0 0 0 0 0 0 13
14 475 101.2 0 0 1 0 0 0 0 0 0 0 0 0 14
15 470 109.6 0 0 0 1 0 0 0 0 0 0 0 0 15
16 461 111.9 0 0 0 0 1 0 0 0 0 0 0 0 16
17 455 108.9 0 0 0 0 0 1 0 0 0 0 0 0 17
18 456 105.6 0 0 0 0 0 0 1 0 0 0 0 0 18
19 517 107.8 0 0 0 0 0 0 0 1 0 0 0 0 19
20 525 97.5 0 0 0 0 0 0 0 0 1 0 0 0 20
21 523 102.4 0 0 0 0 0 0 0 0 0 1 0 0 21
22 519 105.6 0 0 0 0 0 0 0 0 0 0 1 0 22
23 509 99.8 0 0 0 0 0 0 0 0 0 0 0 1 23
24 512 96.2 0 0 0 0 0 0 0 0 0 0 0 0 24
25 519 113.1 0 1 0 0 0 0 0 0 0 0 0 0 25
26 517 107.4 0 0 1 0 0 0 0 0 0 0 0 0 26
27 510 116.8 0 0 0 1 0 0 0 0 0 0 0 0 27
28 509 112.9 0 0 0 0 1 0 0 0 0 0 0 0 28
29 501 105.3 0 0 0 0 0 1 0 0 0 0 0 0 29
30 507 109.3 0 0 0 0 0 0 1 0 0 0 0 0 30
31 569 107.9 0 0 0 0 0 0 0 1 0 0 0 0 31
32 580 101.1 0 0 0 0 0 0 0 0 1 0 0 0 32
33 578 114.7 0 0 0 0 0 0 0 0 0 1 0 0 33
34 565 116.2 0 0 0 0 0 0 0 0 0 0 1 0 34
35 547 108.4 0 0 0 0 0 0 0 0 0 0 0 1 35
36 555 113.4 0 0 0 0 0 0 0 0 0 0 0 0 36
37 562 108.7 0 1 0 0 0 0 0 0 0 0 0 0 37
38 561 112.6 0 0 1 0 0 0 0 0 0 0 0 0 38
39 555 124.2 1 0 0 1 0 0 0 0 0 0 0 0 39
40 544 114.9 1 0 0 0 1 0 0 0 0 0 0 0 40
41 537 110.5 1 0 0 0 0 1 0 0 0 0 0 0 41
42 543 121.5 1 0 0 0 0 0 1 0 0 0 0 0 42
43 594 118.1 1 0 0 0 0 0 0 1 0 0 0 0 43
44 611 111.7 1 0 0 0 0 0 0 0 1 0 0 0 44
45 613 132.7 1 0 0 0 0 0 0 0 0 1 0 0 45
46 611 119.0 1 0 0 0 0 0 0 0 0 0 1 0 46
47 594 116.7 1 0 0 0 0 0 0 0 0 0 0 1 47
48 595 120.1 1 0 0 0 0 0 0 0 0 0 0 0 48
49 591 113.4 1 1 0 0 0 0 0 0 0 0 0 0 49
50 589 106.6 1 0 1 0 0 0 0 0 0 0 0 0 50
51 584 116.3 1 0 0 1 0 0 0 0 0 0 0 0 51
52 573 112.6 1 0 0 0 1 0 0 0 0 0 0 0 52
53 567 111.6 1 0 0 0 0 1 0 0 0 0 0 0 53
54 569 125.1 1 0 0 0 0 0 1 0 0 0 0 0 54
55 621 110.7 1 0 0 0 0 0 0 1 0 0 0 0 55
56 629 109.6 1 0 0 0 0 0 0 0 1 0 0 0 56
57 628 114.2 1 0 0 0 0 0 0 0 0 1 0 0 57
58 612 113.4 1 0 0 0 0 0 0 0 0 0 1 0 58
59 595 116.0 1 0 0 0 0 0 0 0 0 0 0 1 59
60 597 109.6 1 0 0 0 0 0 0 0 0 0 0 0 60
61 593 117.8 1 1 0 0 0 0 0 0 0 0 0 0 61
62 590 115.8 1 0 1 0 0 0 0 0 0 0 0 0 62
63 580 125.3 1 0 0 1 0 0 0 0 0 0 0 0 63
64 574 113.0 1 0 0 0 1 0 0 0 0 0 0 0 64
65 573 120.5 1 0 0 0 0 1 0 0 0 0 0 0 65
66 573 116.6 1 0 0 0 0 0 1 0 0 0 0 0 66
67 620 111.8 1 0 0 0 0 0 0 1 0 0 0 0 67
68 626 115.2 1 0 0 0 0 0 0 0 1 0 0 0 68
69 620 118.6 1 0 0 0 0 0 0 0 0 1 0 0 69
70 588 122.4 1 0 0 0 0 0 0 0 0 0 1 0 70
71 566 116.4 1 0 0 0 0 0 0 0 0 0 0 1 71
72 557 114.5 1 0 0 0 0 0 0 0 0 0 0 0 72
73 561 119.8 1 1 0 0 0 0 0 0 0 0 0 0 73
74 549 115.8 1 0 1 0 0 0 0 0 0 0 0 0 74
75 532 127.8 1 0 0 1 0 0 0 0 0 0 0 0 75
76 526 118.8 1 0 0 0 1 0 0 0 0 0 0 0 76
77 511 119.7 1 0 0 0 0 1 0 0 0 0 0 0 77
78 499 118.6 1 0 0 0 0 0 1 0 0 0 0 0 78
79 555 120.8 1 0 0 0 0 0 0 1 0 0 0 0 79
80 565 115.9 1 0 0 0 0 0 0 0 1 0 0 0 80
81 542 109.7 1 0 0 0 0 0 0 0 0 1 0 0 81
82 527 114.8 1 0 0 0 0 0 0 0 0 0 1 0 82
83 510 116.2 1 0 0 0 0 0 0 0 0 0 0 1 83
84 514 112.2 1 0 0 0 0 0 0 0 0 0 0 0 84
> k <- length(x[1,])
> df <- as.data.frame(x)
> (mylm <- lm(df))
Call:
lm(formula = df)
Coefficients:
(Intercept) X DUM M1 M2 M3
253.1213 2.5547 53.4799 -9.2002 -5.9039 -46.5270
M4 M5 M6 M7 M8 M9
-37.6544 -42.9023 -50.6722 15.2830 36.4342 18.3111
M10 M11 t
3.7422 -6.8667 -0.3257
> (mysum <- summary(mylm))
Call:
lm(formula = df)
Residuals:
Min 1Q Median 3Q Max
-59.550 -21.899 -1.864 22.955 52.440
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 253.1213 65.8521 3.844 0.000266 ***
X 2.5547 0.6702 3.812 0.000296 ***
DUM 53.4799 13.6716 3.912 0.000212 ***
M1 -9.2002 16.7680 -0.549 0.584999
M2 -5.9039 16.3810 -0.360 0.719639
M3 -46.5270 18.2136 -2.555 0.012844 *
M4 -37.6544 16.7599 -2.247 0.027862 *
M5 -42.9023 16.6125 -2.583 0.011932 *
M6 -50.6722 17.0647 -2.969 0.004101 **
M7 15.2830 16.3766 0.933 0.353961
M8 36.4342 16.2152 2.247 0.027847 *
M9 18.3111 16.5529 1.106 0.272476
M10 3.7422 16.6135 0.225 0.822452
M11 -6.8667 16.2485 -0.423 0.673897
t -0.3257 0.3025 -1.077 0.285334
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 30.26 on 69 degrees of freedom
Multiple R-squared: 0.7169, Adjusted R-squared: 0.6594
F-statistic: 12.48 on 14 and 69 DF, p-value: 7.492e-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,] 9.675980e-03 1.935196e-02 0.9903240
[2,] 3.800596e-03 7.601192e-03 0.9961994
[3,] 7.899858e-04 1.579972e-03 0.9992100
[4,] 5.370144e-04 1.074029e-03 0.9994630
[5,] 1.309856e-04 2.619712e-04 0.9998690
[6,] 4.244500e-04 8.489001e-04 0.9995755
[7,] 6.364179e-04 1.272836e-03 0.9993636
[8,] 5.586947e-04 1.117389e-03 0.9994413
[9,] 6.009008e-04 1.201802e-03 0.9993991
[10,] 4.270786e-04 8.541572e-04 0.9995729
[11,] 4.621804e-04 9.243609e-04 0.9995378
[12,] 3.338771e-04 6.677542e-04 0.9996661
[13,] 3.661834e-04 7.323669e-04 0.9996338
[14,] 5.285833e-04 1.057167e-03 0.9994714
[15,] 3.341237e-04 6.682474e-04 0.9996659
[16,] 2.304023e-04 4.608047e-04 0.9997696
[17,] 1.008747e-04 2.017494e-04 0.9998991
[18,] 4.625549e-05 9.251097e-05 0.9999537
[19,] 2.561005e-05 5.122010e-05 0.9999744
[20,] 1.140959e-05 2.281917e-05 0.9999886
[21,] 5.013632e-06 1.002726e-05 0.9999950
[22,] 3.616803e-06 7.233606e-06 0.9999964
[23,] 3.838965e-06 7.677929e-06 0.9999962
[24,] 6.100686e-06 1.220137e-05 0.9999939
[25,] 7.705855e-06 1.541171e-05 0.9999923
[26,] 1.594476e-05 3.188953e-05 0.9999841
[27,] 3.273739e-05 6.547478e-05 0.9999673
[28,] 5.673778e-05 1.134756e-04 0.9999433
[29,] 5.375293e-05 1.075059e-04 0.9999462
[30,] 5.892419e-05 1.178484e-04 0.9999411
[31,] 1.700602e-04 3.401204e-04 0.9998299
[32,] 2.649207e-04 5.298414e-04 0.9997351
[33,] 2.228050e-04 4.456101e-04 0.9997772
[34,] 1.318756e-04 2.637512e-04 0.9998681
[35,] 2.334972e-04 4.669944e-04 0.9997665
[36,] 3.264742e-04 6.529485e-04 0.9996735
[37,] 2.442935e-03 4.885869e-03 0.9975571
[38,] 4.300480e-03 8.600960e-03 0.9956995
[39,] 1.806890e-02 3.613779e-02 0.9819311
[40,] 5.396607e-02 1.079321e-01 0.9460339
[41,] 9.084178e-02 1.816836e-01 0.9091582
[42,] 1.995522e-01 3.991044e-01 0.8004478
[43,] 2.854206e-01 5.708412e-01 0.7145794
[44,] 5.415892e-01 9.168216e-01 0.4584108
[45,] 6.761646e-01 6.476707e-01 0.3238354
[46,] 6.871372e-01 6.257256e-01 0.3128628
[47,] 6.606171e-01 6.787658e-01 0.3393829
[48,] 5.359459e-01 9.281081e-01 0.4640541
[49,] 4.861863e-01 9.723725e-01 0.5138137
> postscript(file="/var/www/html/freestat/rcomp/tmp/1px0q1229363277.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/freestat/rcomp/tmp/228wf1229363277.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/freestat/rcomp/tmp/360dy1229363277.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/freestat/rcomp/tmp/403zr1229363277.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/freestat/rcomp/tmp/54s0j1229363277.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 = 84
Frequency = 1
1 2 3 4 5 6
-34.6151363 -38.7099866 -26.1109604 -22.5416749 -32.3984506 -32.3683850
7 8 9 10 11 12
-17.6630390 -19.3861426 7.4275251 -2.2393664 -12.5752306 -1.4891364
13 14 15 16 17 18
-27.8376196 -26.1879350 -11.6982102 -35.1207813 -27.8832574 -10.3572635
19 20 21 22 23 24
-20.6069118 -7.1195561 -3.1884434 -0.4687081 15.2828688 20.9386156
25 26 27 28 29 30
-5.7090532 3.8819090 13.8169830 14.2332464 31.2221636 35.0992072
31 32 33 34 35 36
35.0463015 42.5923797 24.2980311 22.3606726 35.2215510 23.9073016
37 38 39 40 41 42
52.4400884 38.5064036 -9.6586913 -5.4473140 4.3667209 -9.6387906
43 44 45 46 47 48
-15.5823948 -3.0581769 -36.2569409 11.6363917 11.4466911 -2.7801172
49 50 51 52 53 54
19.8619711 32.2630490 43.4317278 33.3370611 35.4652835 11.0731453
55 56 57 58 59 60
34.2306990 24.2152681 29.9127759 30.8511141 18.1436250 29.9523938
61 62 63 64 65 66
14.5301863 13.6689408 20.3485497 37.2238792 22.6375705 40.6963548
67 68 69 70 71 72
34.3292616 10.8179024 14.5809912 -12.2320640 -7.9695569 -18.6567163
73 74 75 76 77 78
-18.6704367 -23.4223808 -30.1293987 -21.6844166 -33.4100305 -34.5042682
79 80 81 82 83 84
-49.7539165 -48.0616747 -36.7739390 -49.9080400 -59.5499483 -51.8723412
> postscript(file="/var/www/html/freestat/rcomp/tmp/6ubcy1229363277.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 = 84
Frequency = 1
lag(myerror, k = 1) myerror
0 -34.6151363 NA
1 -38.7099866 -34.6151363
2 -26.1109604 -38.7099866
3 -22.5416749 -26.1109604
4 -32.3984506 -22.5416749
5 -32.3683850 -32.3984506
6 -17.6630390 -32.3683850
7 -19.3861426 -17.6630390
8 7.4275251 -19.3861426
9 -2.2393664 7.4275251
10 -12.5752306 -2.2393664
11 -1.4891364 -12.5752306
12 -27.8376196 -1.4891364
13 -26.1879350 -27.8376196
14 -11.6982102 -26.1879350
15 -35.1207813 -11.6982102
16 -27.8832574 -35.1207813
17 -10.3572635 -27.8832574
18 -20.6069118 -10.3572635
19 -7.1195561 -20.6069118
20 -3.1884434 -7.1195561
21 -0.4687081 -3.1884434
22 15.2828688 -0.4687081
23 20.9386156 15.2828688
24 -5.7090532 20.9386156
25 3.8819090 -5.7090532
26 13.8169830 3.8819090
27 14.2332464 13.8169830
28 31.2221636 14.2332464
29 35.0992072 31.2221636
30 35.0463015 35.0992072
31 42.5923797 35.0463015
32 24.2980311 42.5923797
33 22.3606726 24.2980311
34 35.2215510 22.3606726
35 23.9073016 35.2215510
36 52.4400884 23.9073016
37 38.5064036 52.4400884
38 -9.6586913 38.5064036
39 -5.4473140 -9.6586913
40 4.3667209 -5.4473140
41 -9.6387906 4.3667209
42 -15.5823948 -9.6387906
43 -3.0581769 -15.5823948
44 -36.2569409 -3.0581769
45 11.6363917 -36.2569409
46 11.4466911 11.6363917
47 -2.7801172 11.4466911
48 19.8619711 -2.7801172
49 32.2630490 19.8619711
50 43.4317278 32.2630490
51 33.3370611 43.4317278
52 35.4652835 33.3370611
53 11.0731453 35.4652835
54 34.2306990 11.0731453
55 24.2152681 34.2306990
56 29.9127759 24.2152681
57 30.8511141 29.9127759
58 18.1436250 30.8511141
59 29.9523938 18.1436250
60 14.5301863 29.9523938
61 13.6689408 14.5301863
62 20.3485497 13.6689408
63 37.2238792 20.3485497
64 22.6375705 37.2238792
65 40.6963548 22.6375705
66 34.3292616 40.6963548
67 10.8179024 34.3292616
68 14.5809912 10.8179024
69 -12.2320640 14.5809912
70 -7.9695569 -12.2320640
71 -18.6567163 -7.9695569
72 -18.6704367 -18.6567163
73 -23.4223808 -18.6704367
74 -30.1293987 -23.4223808
75 -21.6844166 -30.1293987
76 -33.4100305 -21.6844166
77 -34.5042682 -33.4100305
78 -49.7539165 -34.5042682
79 -48.0616747 -49.7539165
80 -36.7739390 -48.0616747
81 -49.9080400 -36.7739390
82 -59.5499483 -49.9080400
83 -51.8723412 -59.5499483
84 NA -51.8723412
> dum1 <- dum[2:length(myerror),]
> dum1
lag(myerror, k = 1) myerror
[1,] -38.7099866 -34.6151363
[2,] -26.1109604 -38.7099866
[3,] -22.5416749 -26.1109604
[4,] -32.3984506 -22.5416749
[5,] -32.3683850 -32.3984506
[6,] -17.6630390 -32.3683850
[7,] -19.3861426 -17.6630390
[8,] 7.4275251 -19.3861426
[9,] -2.2393664 7.4275251
[10,] -12.5752306 -2.2393664
[11,] -1.4891364 -12.5752306
[12,] -27.8376196 -1.4891364
[13,] -26.1879350 -27.8376196
[14,] -11.6982102 -26.1879350
[15,] -35.1207813 -11.6982102
[16,] -27.8832574 -35.1207813
[17,] -10.3572635 -27.8832574
[18,] -20.6069118 -10.3572635
[19,] -7.1195561 -20.6069118
[20,] -3.1884434 -7.1195561
[21,] -0.4687081 -3.1884434
[22,] 15.2828688 -0.4687081
[23,] 20.9386156 15.2828688
[24,] -5.7090532 20.9386156
[25,] 3.8819090 -5.7090532
[26,] 13.8169830 3.8819090
[27,] 14.2332464 13.8169830
[28,] 31.2221636 14.2332464
[29,] 35.0992072 31.2221636
[30,] 35.0463015 35.0992072
[31,] 42.5923797 35.0463015
[32,] 24.2980311 42.5923797
[33,] 22.3606726 24.2980311
[34,] 35.2215510 22.3606726
[35,] 23.9073016 35.2215510
[36,] 52.4400884 23.9073016
[37,] 38.5064036 52.4400884
[38,] -9.6586913 38.5064036
[39,] -5.4473140 -9.6586913
[40,] 4.3667209 -5.4473140
[41,] -9.6387906 4.3667209
[42,] -15.5823948 -9.6387906
[43,] -3.0581769 -15.5823948
[44,] -36.2569409 -3.0581769
[45,] 11.6363917 -36.2569409
[46,] 11.4466911 11.6363917
[47,] -2.7801172 11.4466911
[48,] 19.8619711 -2.7801172
[49,] 32.2630490 19.8619711
[50,] 43.4317278 32.2630490
[51,] 33.3370611 43.4317278
[52,] 35.4652835 33.3370611
[53,] 11.0731453 35.4652835
[54,] 34.2306990 11.0731453
[55,] 24.2152681 34.2306990
[56,] 29.9127759 24.2152681
[57,] 30.8511141 29.9127759
[58,] 18.1436250 30.8511141
[59,] 29.9523938 18.1436250
[60,] 14.5301863 29.9523938
[61,] 13.6689408 14.5301863
[62,] 20.3485497 13.6689408
[63,] 37.2238792 20.3485497
[64,] 22.6375705 37.2238792
[65,] 40.6963548 22.6375705
[66,] 34.3292616 40.6963548
[67,] 10.8179024 34.3292616
[68,] 14.5809912 10.8179024
[69,] -12.2320640 14.5809912
[70,] -7.9695569 -12.2320640
[71,] -18.6567163 -7.9695569
[72,] -18.6704367 -18.6567163
[73,] -23.4223808 -18.6704367
[74,] -30.1293987 -23.4223808
[75,] -21.6844166 -30.1293987
[76,] -33.4100305 -21.6844166
[77,] -34.5042682 -33.4100305
[78,] -49.7539165 -34.5042682
[79,] -48.0616747 -49.7539165
[80,] -36.7739390 -48.0616747
[81,] -49.9080400 -36.7739390
[82,] -59.5499483 -49.9080400
[83,] -51.8723412 -59.5499483
> z <- as.data.frame(dum1)
> z
lag(myerror, k = 1) myerror
1 -38.7099866 -34.6151363
2 -26.1109604 -38.7099866
3 -22.5416749 -26.1109604
4 -32.3984506 -22.5416749
5 -32.3683850 -32.3984506
6 -17.6630390 -32.3683850
7 -19.3861426 -17.6630390
8 7.4275251 -19.3861426
9 -2.2393664 7.4275251
10 -12.5752306 -2.2393664
11 -1.4891364 -12.5752306
12 -27.8376196 -1.4891364
13 -26.1879350 -27.8376196
14 -11.6982102 -26.1879350
15 -35.1207813 -11.6982102
16 -27.8832574 -35.1207813
17 -10.3572635 -27.8832574
18 -20.6069118 -10.3572635
19 -7.1195561 -20.6069118
20 -3.1884434 -7.1195561
21 -0.4687081 -3.1884434
22 15.2828688 -0.4687081
23 20.9386156 15.2828688
24 -5.7090532 20.9386156
25 3.8819090 -5.7090532
26 13.8169830 3.8819090
27 14.2332464 13.8169830
28 31.2221636 14.2332464
29 35.0992072 31.2221636
30 35.0463015 35.0992072
31 42.5923797 35.0463015
32 24.2980311 42.5923797
33 22.3606726 24.2980311
34 35.2215510 22.3606726
35 23.9073016 35.2215510
36 52.4400884 23.9073016
37 38.5064036 52.4400884
38 -9.6586913 38.5064036
39 -5.4473140 -9.6586913
40 4.3667209 -5.4473140
41 -9.6387906 4.3667209
42 -15.5823948 -9.6387906
43 -3.0581769 -15.5823948
44 -36.2569409 -3.0581769
45 11.6363917 -36.2569409
46 11.4466911 11.6363917
47 -2.7801172 11.4466911
48 19.8619711 -2.7801172
49 32.2630490 19.8619711
50 43.4317278 32.2630490
51 33.3370611 43.4317278
52 35.4652835 33.3370611
53 11.0731453 35.4652835
54 34.2306990 11.0731453
55 24.2152681 34.2306990
56 29.9127759 24.2152681
57 30.8511141 29.9127759
58 18.1436250 30.8511141
59 29.9523938 18.1436250
60 14.5301863 29.9523938
61 13.6689408 14.5301863
62 20.3485497 13.6689408
63 37.2238792 20.3485497
64 22.6375705 37.2238792
65 40.6963548 22.6375705
66 34.3292616 40.6963548
67 10.8179024 34.3292616
68 14.5809912 10.8179024
69 -12.2320640 14.5809912
70 -7.9695569 -12.2320640
71 -18.6567163 -7.9695569
72 -18.6704367 -18.6567163
73 -23.4223808 -18.6704367
74 -30.1293987 -23.4223808
75 -21.6844166 -30.1293987
76 -33.4100305 -21.6844166
77 -34.5042682 -33.4100305
78 -49.7539165 -34.5042682
79 -48.0616747 -49.7539165
80 -36.7739390 -48.0616747
81 -49.9080400 -36.7739390
82 -59.5499483 -49.9080400
83 -51.8723412 -59.5499483
> 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/freestat/rcomp/tmp/7o1t21229363277.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/freestat/rcomp/tmp/81tv11229363277.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/freestat/rcomp/tmp/9ncav1229363277.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/freestat/rcomp/tmp/10qir31229363277.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/freestat/rcomp/createtable file can be downloaded at http://www.wessa.net/cretab
> load(file="/var/www/html/freestat/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/freestat/rcomp/tmp/11pqoz1229363277.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/freestat/rcomp/tmp/12jkhk1229363277.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/freestat/rcomp/tmp/13lg2l1229363277.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/freestat/rcomp/tmp/14aijn1229363277.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/freestat/rcomp/tmp/15isw81229363277.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/freestat/rcomp/tmp/166zmf1229363278.tab")
+ }
>
> system("convert tmp/1px0q1229363277.ps tmp/1px0q1229363277.png")
> system("convert tmp/228wf1229363277.ps tmp/228wf1229363277.png")
> system("convert tmp/360dy1229363277.ps tmp/360dy1229363277.png")
> system("convert tmp/403zr1229363277.ps tmp/403zr1229363277.png")
> system("convert tmp/54s0j1229363277.ps tmp/54s0j1229363277.png")
> system("convert tmp/6ubcy1229363277.ps tmp/6ubcy1229363277.png")
> system("convert tmp/7o1t21229363277.ps tmp/7o1t21229363277.png")
> system("convert tmp/81tv11229363277.ps tmp/81tv11229363277.png")
> system("convert tmp/9ncav1229363277.ps tmp/9ncav1229363277.png")
> system("convert tmp/10qir31229363277.ps tmp/10qir31229363277.png")
>
>
> proc.time()
user system elapsed
4.046 2.534 4.971