R version 2.15.2 (2012-10-26) -- "Trick or Treat"
Copyright (C) 2012 The R Foundation for Statistical Computing
ISBN 3-900051-07-0
Platform: i686-pc-linux-gnu (32-bit)
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(1.2999,1.3074,1.3242,1.3516,1.3511,1.3419,1.3716,1.3622,1.3896,1.4227,1.4684,1.457,1.4718,1.4748,1.5527,1.5751,1.5557,1.5553,1.577,1.4975,1.437,1.3322,1.2732,1.3449,1.3239,1.2785,1.305,1.319,1.365,1.4016,1.4088,1.4268,1.4562,1.4816,1.4914,1.4614,1.4272,1.3686,1.3569,1.3406,1.2565,1.2209,1.277,1.2894,1.3067,1.3898,1.3661,1.322,1.336,1.3649,1.3999,1.4442,1.4349,1.4388,1.4264,1.4343,1.377,1.3706,1.3556,1.3179,1.2905,1.3224,1.3201,1.3162,1.2789,1.2526,1.2288,1.24,1.2856),dim=c(1,69),dimnames=list(c('Exchange_rate'),1:69))
> y <- array(NA,dim=c(1,69),dimnames=list(c('Exchange_rate'),1:69))
> 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'
> library(lattice)
> library(lmtest)
Loading required package: zoo
Attaching package: 'zoo'
The following object(s) are masked from 'package:base':
as.Date, 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
Exchange_rate M1 M2 M3 M4 M5 M6 M7 M8 M9 M10 M11 t
1 1.2999 1 0 0 0 0 0 0 0 0 0 0 1
2 1.3074 0 1 0 0 0 0 0 0 0 0 0 2
3 1.3242 0 0 1 0 0 0 0 0 0 0 0 3
4 1.3516 0 0 0 1 0 0 0 0 0 0 0 4
5 1.3511 0 0 0 0 1 0 0 0 0 0 0 5
6 1.3419 0 0 0 0 0 1 0 0 0 0 0 6
7 1.3716 0 0 0 0 0 0 1 0 0 0 0 7
8 1.3622 0 0 0 0 0 0 0 1 0 0 0 8
9 1.3896 0 0 0 0 0 0 0 0 1 0 0 9
10 1.4227 0 0 0 0 0 0 0 0 0 1 0 10
11 1.4684 0 0 0 0 0 0 0 0 0 0 1 11
12 1.4570 0 0 0 0 0 0 0 0 0 0 0 12
13 1.4718 1 0 0 0 0 0 0 0 0 0 0 13
14 1.4748 0 1 0 0 0 0 0 0 0 0 0 14
15 1.5527 0 0 1 0 0 0 0 0 0 0 0 15
16 1.5751 0 0 0 1 0 0 0 0 0 0 0 16
17 1.5557 0 0 0 0 1 0 0 0 0 0 0 17
18 1.5553 0 0 0 0 0 1 0 0 0 0 0 18
19 1.5770 0 0 0 0 0 0 1 0 0 0 0 19
20 1.4975 0 0 0 0 0 0 0 1 0 0 0 20
21 1.4370 0 0 0 0 0 0 0 0 1 0 0 21
22 1.3322 0 0 0 0 0 0 0 0 0 1 0 22
23 1.2732 0 0 0 0 0 0 0 0 0 0 1 23
24 1.3449 0 0 0 0 0 0 0 0 0 0 0 24
25 1.3239 1 0 0 0 0 0 0 0 0 0 0 25
26 1.2785 0 1 0 0 0 0 0 0 0 0 0 26
27 1.3050 0 0 1 0 0 0 0 0 0 0 0 27
28 1.3190 0 0 0 1 0 0 0 0 0 0 0 28
29 1.3650 0 0 0 0 1 0 0 0 0 0 0 29
30 1.4016 0 0 0 0 0 1 0 0 0 0 0 30
31 1.4088 0 0 0 0 0 0 1 0 0 0 0 31
32 1.4268 0 0 0 0 0 0 0 1 0 0 0 32
33 1.4562 0 0 0 0 0 0 0 0 1 0 0 33
34 1.4816 0 0 0 0 0 0 0 0 0 1 0 34
35 1.4914 0 0 0 0 0 0 0 0 0 0 1 35
36 1.4614 0 0 0 0 0 0 0 0 0 0 0 36
37 1.4272 1 0 0 0 0 0 0 0 0 0 0 37
38 1.3686 0 1 0 0 0 0 0 0 0 0 0 38
39 1.3569 0 0 1 0 0 0 0 0 0 0 0 39
40 1.3406 0 0 0 1 0 0 0 0 0 0 0 40
41 1.2565 0 0 0 0 1 0 0 0 0 0 0 41
42 1.2209 0 0 0 0 0 1 0 0 0 0 0 42
43 1.2770 0 0 0 0 0 0 1 0 0 0 0 43
44 1.2894 0 0 0 0 0 0 0 1 0 0 0 44
45 1.3067 0 0 0 0 0 0 0 0 1 0 0 45
46 1.3898 0 0 0 0 0 0 0 0 0 1 0 46
47 1.3661 0 0 0 0 0 0 0 0 0 0 1 47
48 1.3220 0 0 0 0 0 0 0 0 0 0 0 48
49 1.3360 1 0 0 0 0 0 0 0 0 0 0 49
50 1.3649 0 1 0 0 0 0 0 0 0 0 0 50
51 1.3999 0 0 1 0 0 0 0 0 0 0 0 51
52 1.4442 0 0 0 1 0 0 0 0 0 0 0 52
53 1.4349 0 0 0 0 1 0 0 0 0 0 0 53
54 1.4388 0 0 0 0 0 1 0 0 0 0 0 54
55 1.4264 0 0 0 0 0 0 1 0 0 0 0 55
56 1.4343 0 0 0 0 0 0 0 1 0 0 0 56
57 1.3770 0 0 0 0 0 0 0 0 1 0 0 57
58 1.3706 0 0 0 0 0 0 0 0 0 1 0 58
59 1.3556 0 0 0 0 0 0 0 0 0 0 1 59
60 1.3179 0 0 0 0 0 0 0 0 0 0 0 60
61 1.2905 1 0 0 0 0 0 0 0 0 0 0 61
62 1.3224 0 1 0 0 0 0 0 0 0 0 0 62
63 1.3201 0 0 1 0 0 0 0 0 0 0 0 63
64 1.3162 0 0 0 1 0 0 0 0 0 0 0 64
65 1.2789 0 0 0 0 1 0 0 0 0 0 0 65
66 1.2526 0 0 0 0 0 1 0 0 0 0 0 66
67 1.2288 0 0 0 0 0 0 1 0 0 0 0 67
68 1.2400 0 0 0 0 0 0 0 1 0 0 0 68
69 1.2856 0 0 0 0 0 0 0 0 1 0 0 69
> k <- length(x[1,])
> df <- as.data.frame(x)
> (mylm <- lm(df))
Call:
lm(formula = df)
Coefficients:
(Intercept) M1 M2 M3 M4 M5
1.4394024 -0.0305848 -0.0344025 -0.0090702 0.0072121 -0.0085890
M6 M7 M8 M9 M10 M11
-0.0121233 0.0025923 -0.0023421 -0.0003931 0.0154754 0.0086677
t
-0.0016323
> (mysum <- summary(mylm))
Call:
lm(formula = df)
Residuals:
Min 1Q Median 3Q Max
-0.137823 -0.066948 -0.007398 0.071056 0.166019
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 1.4394024 0.0432260 33.299 < 2e-16 ***
M1 -0.0305848 0.0526263 -0.581 0.56346
M2 -0.0344025 0.0526024 -0.654 0.51578
M3 -0.0090702 0.0525839 -0.172 0.86367
M4 0.0072121 0.0525706 0.137 0.89137
M5 -0.0085890 0.0525626 -0.163 0.87079
M6 -0.0121233 0.0525600 -0.231 0.81842
M7 0.0025923 0.0525626 0.049 0.96084
M8 -0.0023421 0.0525706 -0.045 0.96462
M9 -0.0003931 0.0525839 -0.007 0.99406
M10 0.0154754 0.0549073 0.282 0.77910
M11 0.0086677 0.0548997 0.158 0.87512
t -0.0016323 0.0005282 -3.090 0.00311 **
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 0.0868 on 56 degrees of freedom
Multiple R-squared: 0.1643, Adjusted R-squared: -0.01476
F-statistic: 0.9176 on 12 and 56 DF, p-value: 0.5359
> 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.0270476714 0.054095343 0.972952329
[2,] 0.0068142079 0.013628416 0.993185792
[3,] 0.0019109288 0.003821858 0.998089071
[4,] 0.0005502111 0.001100422 0.999449789
[5,] 0.0015322612 0.003064522 0.998467739
[6,] 0.0333109913 0.066621983 0.966689009
[7,] 0.4541857363 0.908371473 0.545814264
[8,] 0.9201062961 0.159787408 0.079893704
[9,] 0.9529829138 0.094034172 0.047017086
[10,] 0.9699639451 0.060072110 0.030036055
[11,] 0.9857607732 0.028478454 0.014239227
[12,] 0.9918212153 0.016357569 0.008178785
[13,] 0.9947702549 0.010459490 0.005229745
[14,] 0.9922186400 0.015562720 0.007781360
[15,] 0.9868268464 0.026346307 0.013173154
[16,] 0.9790807191 0.041838562 0.020919281
[17,] 0.9665030309 0.066993938 0.033496969
[18,] 0.9544262337 0.091147533 0.045573766
[19,] 0.9424285424 0.115142915 0.057571458
[20,] 0.9375910514 0.124817897 0.062408949
[21,] 0.9311131590 0.137773682 0.068886841
[22,] 0.9167336536 0.166532693 0.083266346
[23,] 0.8766480333 0.246703933 0.123351967
[24,] 0.8312458043 0.337508391 0.168754196
[25,] 0.7973116530 0.405376694 0.202688347
[26,] 0.8478407035 0.304318593 0.152159297
[27,] 0.9349851092 0.130029782 0.065014891
[28,] 0.9521513222 0.095697356 0.047848678
[29,] 0.9726333650 0.054733270 0.027366635
[30,] 0.9900541142 0.019891772 0.009945886
[31,] 0.9863047886 0.027390423 0.013695211
[32,] 0.9851163224 0.029767355 0.014883678
[33,] 0.9895898445 0.020820311 0.010410156
[34,] 0.9880849474 0.023830105 0.011915053
[35,] 0.9921978739 0.015604252 0.007802126
[36,] 0.9920159880 0.015968024 0.007984012
[37,] 0.9793095658 0.041380868 0.020690434
[38,] 0.9340355535 0.131928893 0.065964446
> postscript(file="/var/wessaorg/rcomp/tmp/1fi4b1353793750.ps",horizontal=F,onefile=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/wessaorg/rcomp/tmp/25g8p1353793750.ps",horizontal=F,onefile=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/wessaorg/rcomp/tmp/3gr361353793750.ps",horizontal=F,onefile=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/wessaorg/rcomp/tmp/4wqch1353793750.ps",horizontal=F,onefile=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/wessaorg/rcomp/tmp/59ym31353793750.ps",horizontal=F,onefile=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 = 69
Frequency = 1
1 2 3 4 5 6
-0.107285333 -0.094335333 -0.101235333 -0.088485333 -0.071552000 -0.075585333
7 8 9 10 11 12
-0.058968667 -0.061802000 -0.034718667 -0.015854933 0.038285067 0.037185067
13 14 15 16 17 18
0.084202133 0.092652133 0.146852133 0.154602133 0.152635467 0.157402133
19 20 21 22 23 24
0.166018800 0.093085467 0.032268800 -0.086767467 -0.137327467 -0.055327467
25 26 27 28 29 30
-0.044110400 -0.084060400 -0.081260400 -0.081910400 -0.018477067 0.023289600
31 32 33 34 35 36
0.017406267 0.041972933 0.071056267 0.082220000 0.100460000 0.080760000
37 38 39 40 41 42
0.078777067 0.025627067 -0.009772933 -0.040722933 -0.107389600 -0.137822933
43 44 45 46 47 48
-0.094806267 -0.075839600 -0.058856267 0.010007467 -0.005252533 -0.039052533
49 50 51 52 53 54
0.007164533 0.041514533 0.052814533 0.082464533 0.090597867 0.099664533
55 56 57 58 59 60
0.074181200 0.088647867 0.031031200 0.010394933 0.003834933 -0.023565067
61 62 63 64 65 66
-0.018748000 0.018602000 -0.007398000 -0.025948000 -0.045814667 -0.066948000
67 68 69
-0.103831333 -0.086064667 -0.040781333
> postscript(file="/var/wessaorg/rcomp/tmp/6q3y81353793750.ps",horizontal=F,onefile=F,pagecentre=F,paper="special",width=8.3333333333333,height=5.5555555555556)
> dum <- cbind(lag(myerror,k=1),myerror)
> dum
Time Series:
Start = 0
End = 69
Frequency = 1
lag(myerror, k = 1) myerror
0 -0.107285333 NA
1 -0.094335333 -0.107285333
2 -0.101235333 -0.094335333
3 -0.088485333 -0.101235333
4 -0.071552000 -0.088485333
5 -0.075585333 -0.071552000
6 -0.058968667 -0.075585333
7 -0.061802000 -0.058968667
8 -0.034718667 -0.061802000
9 -0.015854933 -0.034718667
10 0.038285067 -0.015854933
11 0.037185067 0.038285067
12 0.084202133 0.037185067
13 0.092652133 0.084202133
14 0.146852133 0.092652133
15 0.154602133 0.146852133
16 0.152635467 0.154602133
17 0.157402133 0.152635467
18 0.166018800 0.157402133
19 0.093085467 0.166018800
20 0.032268800 0.093085467
21 -0.086767467 0.032268800
22 -0.137327467 -0.086767467
23 -0.055327467 -0.137327467
24 -0.044110400 -0.055327467
25 -0.084060400 -0.044110400
26 -0.081260400 -0.084060400
27 -0.081910400 -0.081260400
28 -0.018477067 -0.081910400
29 0.023289600 -0.018477067
30 0.017406267 0.023289600
31 0.041972933 0.017406267
32 0.071056267 0.041972933
33 0.082220000 0.071056267
34 0.100460000 0.082220000
35 0.080760000 0.100460000
36 0.078777067 0.080760000
37 0.025627067 0.078777067
38 -0.009772933 0.025627067
39 -0.040722933 -0.009772933
40 -0.107389600 -0.040722933
41 -0.137822933 -0.107389600
42 -0.094806267 -0.137822933
43 -0.075839600 -0.094806267
44 -0.058856267 -0.075839600
45 0.010007467 -0.058856267
46 -0.005252533 0.010007467
47 -0.039052533 -0.005252533
48 0.007164533 -0.039052533
49 0.041514533 0.007164533
50 0.052814533 0.041514533
51 0.082464533 0.052814533
52 0.090597867 0.082464533
53 0.099664533 0.090597867
54 0.074181200 0.099664533
55 0.088647867 0.074181200
56 0.031031200 0.088647867
57 0.010394933 0.031031200
58 0.003834933 0.010394933
59 -0.023565067 0.003834933
60 -0.018748000 -0.023565067
61 0.018602000 -0.018748000
62 -0.007398000 0.018602000
63 -0.025948000 -0.007398000
64 -0.045814667 -0.025948000
65 -0.066948000 -0.045814667
66 -0.103831333 -0.066948000
67 -0.086064667 -0.103831333
68 -0.040781333 -0.086064667
69 NA -0.040781333
> dum1 <- dum[2:length(myerror),]
> dum1
lag(myerror, k = 1) myerror
[1,] -0.094335333 -0.107285333
[2,] -0.101235333 -0.094335333
[3,] -0.088485333 -0.101235333
[4,] -0.071552000 -0.088485333
[5,] -0.075585333 -0.071552000
[6,] -0.058968667 -0.075585333
[7,] -0.061802000 -0.058968667
[8,] -0.034718667 -0.061802000
[9,] -0.015854933 -0.034718667
[10,] 0.038285067 -0.015854933
[11,] 0.037185067 0.038285067
[12,] 0.084202133 0.037185067
[13,] 0.092652133 0.084202133
[14,] 0.146852133 0.092652133
[15,] 0.154602133 0.146852133
[16,] 0.152635467 0.154602133
[17,] 0.157402133 0.152635467
[18,] 0.166018800 0.157402133
[19,] 0.093085467 0.166018800
[20,] 0.032268800 0.093085467
[21,] -0.086767467 0.032268800
[22,] -0.137327467 -0.086767467
[23,] -0.055327467 -0.137327467
[24,] -0.044110400 -0.055327467
[25,] -0.084060400 -0.044110400
[26,] -0.081260400 -0.084060400
[27,] -0.081910400 -0.081260400
[28,] -0.018477067 -0.081910400
[29,] 0.023289600 -0.018477067
[30,] 0.017406267 0.023289600
[31,] 0.041972933 0.017406267
[32,] 0.071056267 0.041972933
[33,] 0.082220000 0.071056267
[34,] 0.100460000 0.082220000
[35,] 0.080760000 0.100460000
[36,] 0.078777067 0.080760000
[37,] 0.025627067 0.078777067
[38,] -0.009772933 0.025627067
[39,] -0.040722933 -0.009772933
[40,] -0.107389600 -0.040722933
[41,] -0.137822933 -0.107389600
[42,] -0.094806267 -0.137822933
[43,] -0.075839600 -0.094806267
[44,] -0.058856267 -0.075839600
[45,] 0.010007467 -0.058856267
[46,] -0.005252533 0.010007467
[47,] -0.039052533 -0.005252533
[48,] 0.007164533 -0.039052533
[49,] 0.041514533 0.007164533
[50,] 0.052814533 0.041514533
[51,] 0.082464533 0.052814533
[52,] 0.090597867 0.082464533
[53,] 0.099664533 0.090597867
[54,] 0.074181200 0.099664533
[55,] 0.088647867 0.074181200
[56,] 0.031031200 0.088647867
[57,] 0.010394933 0.031031200
[58,] 0.003834933 0.010394933
[59,] -0.023565067 0.003834933
[60,] -0.018748000 -0.023565067
[61,] 0.018602000 -0.018748000
[62,] -0.007398000 0.018602000
[63,] -0.025948000 -0.007398000
[64,] -0.045814667 -0.025948000
[65,] -0.066948000 -0.045814667
[66,] -0.103831333 -0.066948000
[67,] -0.086064667 -0.103831333
[68,] -0.040781333 -0.086064667
> z <- as.data.frame(dum1)
> z
lag(myerror, k = 1) myerror
1 -0.094335333 -0.107285333
2 -0.101235333 -0.094335333
3 -0.088485333 -0.101235333
4 -0.071552000 -0.088485333
5 -0.075585333 -0.071552000
6 -0.058968667 -0.075585333
7 -0.061802000 -0.058968667
8 -0.034718667 -0.061802000
9 -0.015854933 -0.034718667
10 0.038285067 -0.015854933
11 0.037185067 0.038285067
12 0.084202133 0.037185067
13 0.092652133 0.084202133
14 0.146852133 0.092652133
15 0.154602133 0.146852133
16 0.152635467 0.154602133
17 0.157402133 0.152635467
18 0.166018800 0.157402133
19 0.093085467 0.166018800
20 0.032268800 0.093085467
21 -0.086767467 0.032268800
22 -0.137327467 -0.086767467
23 -0.055327467 -0.137327467
24 -0.044110400 -0.055327467
25 -0.084060400 -0.044110400
26 -0.081260400 -0.084060400
27 -0.081910400 -0.081260400
28 -0.018477067 -0.081910400
29 0.023289600 -0.018477067
30 0.017406267 0.023289600
31 0.041972933 0.017406267
32 0.071056267 0.041972933
33 0.082220000 0.071056267
34 0.100460000 0.082220000
35 0.080760000 0.100460000
36 0.078777067 0.080760000
37 0.025627067 0.078777067
38 -0.009772933 0.025627067
39 -0.040722933 -0.009772933
40 -0.107389600 -0.040722933
41 -0.137822933 -0.107389600
42 -0.094806267 -0.137822933
43 -0.075839600 -0.094806267
44 -0.058856267 -0.075839600
45 0.010007467 -0.058856267
46 -0.005252533 0.010007467
47 -0.039052533 -0.005252533
48 0.007164533 -0.039052533
49 0.041514533 0.007164533
50 0.052814533 0.041514533
51 0.082464533 0.052814533
52 0.090597867 0.082464533
53 0.099664533 0.090597867
54 0.074181200 0.099664533
55 0.088647867 0.074181200
56 0.031031200 0.088647867
57 0.010394933 0.031031200
58 0.003834933 0.010394933
59 -0.023565067 0.003834933
60 -0.018748000 -0.023565067
61 0.018602000 -0.018748000
62 -0.007398000 0.018602000
63 -0.025948000 -0.007398000
64 -0.045814667 -0.025948000
65 -0.066948000 -0.045814667
66 -0.103831333 -0.066948000
67 -0.086064667 -0.103831333
68 -0.040781333 -0.086064667
> 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/wessaorg/rcomp/tmp/7m8h81353793750.ps",horizontal=F,onefile=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/wessaorg/rcomp/tmp/8w9yp1353793750.ps",horizontal=F,onefile=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/wessaorg/rcomp/tmp/9jrxm1353793750.ps",horizontal=F,onefile=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/wessaorg/rcomp/tmp/10r25u1353793750.ps",horizontal=F,onefile=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/wessaorg/rcomp/createtable file can be downloaded at http://www.wessa.net/cretab
> load(file="/var/wessaorg/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/wessaorg/rcomp/tmp/112dq01353793750.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/wessaorg/rcomp/tmp/12c98r1353793750.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/wessaorg/rcomp/tmp/13vf3x1353793750.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/wessaorg/rcomp/tmp/145huo1353793750.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/wessaorg/rcomp/tmp/15zpkc1353793750.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/wessaorg/rcomp/tmp/16ykkk1353793750.tab")
+ }
>
> try(system("convert tmp/1fi4b1353793750.ps tmp/1fi4b1353793750.png",intern=TRUE))
character(0)
> try(system("convert tmp/25g8p1353793750.ps tmp/25g8p1353793750.png",intern=TRUE))
character(0)
> try(system("convert tmp/3gr361353793750.ps tmp/3gr361353793750.png",intern=TRUE))
character(0)
> try(system("convert tmp/4wqch1353793750.ps tmp/4wqch1353793750.png",intern=TRUE))
character(0)
> try(system("convert tmp/59ym31353793750.ps tmp/59ym31353793750.png",intern=TRUE))
character(0)
> try(system("convert tmp/6q3y81353793750.ps tmp/6q3y81353793750.png",intern=TRUE))
character(0)
> try(system("convert tmp/7m8h81353793750.ps tmp/7m8h81353793750.png",intern=TRUE))
character(0)
> try(system("convert tmp/8w9yp1353793750.ps tmp/8w9yp1353793750.png",intern=TRUE))
character(0)
> try(system("convert tmp/9jrxm1353793750.ps tmp/9jrxm1353793750.png",intern=TRUE))
character(0)
> try(system("convert tmp/10r25u1353793750.ps tmp/10r25u1353793750.png",intern=TRUE))
character(0)
>
>
> proc.time()
user system elapsed
5.983 0.841 6.864