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
+ ,591
+ ,13.119
+ ,0.69867
+ ,135.63
+ ,2
+ ,589
+ ,13.014
+ ,0.68968
+ ,136.55
+ ,3
+ ,584
+ ,13.201
+ ,0.69233
+ ,138.83
+ ,4
+ ,573
+ ,12.938
+ ,0.68293
+ ,138.84
+ ,5
+ ,567
+ ,12.694
+ ,0.68399
+ ,135.37
+ ,6
+ ,569
+ ,12.165
+ ,0.66895
+ ,132.22
+ ,7
+ ,621
+ ,12.037
+ ,0.68756
+ ,134.75
+ ,8
+ ,629
+ ,12.292
+ ,0.68527
+ ,135.98
+ ,9
+ ,628
+ ,12.256
+ ,0.6776
+ ,136.06
+ ,10
+ ,612
+ ,12.015
+ ,0.68137
+ ,138.05
+ ,11
+ ,595
+ ,11.786
+ ,0.67933
+ ,139.59
+ ,12
+ ,597
+ ,11.856
+ ,0.67922
+ ,140.58
+ ,13
+ ,593
+ ,12.103
+ ,0.68598
+ ,139.81
+ ,14
+ ,590
+ ,11.938
+ ,0.68297
+ ,140.77
+ ,15
+ ,580
+ ,1.202
+ ,0.68935
+ ,140.96
+ ,16
+ ,574
+ ,12.271
+ ,0.69463
+ ,143.59
+ ,17
+ ,573
+ ,1.277
+ ,0.6833
+ ,142.7
+ ,18
+ ,573
+ ,1.265
+ ,0.68666
+ ,145.11
+ ,19
+ ,620
+ ,12.684
+ ,0.68782
+ ,146.7
+ ,20
+ ,626
+ ,12.811
+ ,0.67669
+ ,148.53
+ ,21
+ ,620
+ ,12.727
+ ,0.67511
+ ,148.99
+ ,22
+ ,588
+ ,12.611
+ ,0.67254
+ ,149.65
+ ,23
+ ,566
+ ,12.881
+ ,0.67397
+ ,151.11
+ ,24
+ ,557
+ ,13.213
+ ,0.67286
+ ,154.82
+ ,25
+ ,561
+ ,12.999
+ ,0.66341
+ ,156.56
+ ,26
+ ,549
+ ,13.074
+ ,0.668
+ ,157.6
+ ,27
+ ,532
+ ,13.242
+ ,0.68021
+ ,155.24
+ ,28
+ ,526
+ ,13.516
+ ,0.67934
+ ,160.68
+ ,29
+ ,511
+ ,13.511
+ ,0.68136
+ ,163.22
+ ,30
+ ,499
+ ,13.419
+ ,0.67562
+ ,164.55
+ ,31
+ ,555
+ ,13.716
+ ,0.6744
+ ,166.76
+ ,32
+ ,565
+ ,13.622
+ ,0.67766
+ ,159.05
+ ,33
+ ,542
+ ,13.896
+ ,0.68887
+ ,159.82
+ ,34
+ ,527
+ ,14.227
+ ,0.69614
+ ,164.95
+ ,35
+ ,510
+ ,14.684
+ ,0.70896
+ ,162.89
+ ,36
+ ,514
+ ,1.457
+ ,0.72064
+ ,163.55
+ ,37
+ ,517
+ ,14.718
+ ,0.74725
+ ,158.68
+ ,38
+ ,508
+ ,14.748
+ ,0.75094
+ ,157.97
+ ,39
+ ,493
+ ,15.527
+ ,0.77494
+ ,156.59
+ ,40
+ ,490
+ ,1.575
+ ,0.79487
+ ,161.56
+ ,41
+ ,469
+ ,15.557
+ ,0.79209
+ ,162.31
+ ,42
+ ,478
+ ,15.553
+ ,0.79152
+ ,166.26
+ ,43
+ ,528
+ ,1.577
+ ,0.79308
+ ,168.45
+ ,44
+ ,534
+ ,14.975
+ ,0.79279
+ ,163.63
+ ,45
+ ,518
+ ,14.369
+ ,0.79924
+ ,153.2
+ ,46
+ ,506
+ ,13.322
+ ,0.78668
+ ,133.52
+ ,47
+ ,502
+ ,12.732
+ ,0.83063
+ ,123.28
+ ,48
+ ,516
+ ,13.449
+ ,0.90448
+ ,122.51
+ ,49
+ ,528
+ ,13.239
+ ,0.91819
+ ,119.73
+ ,50
+ ,533
+ ,12.785
+ ,0.88691
+ ,118.3
+ ,51
+ ,536
+ ,1.305
+ ,0.91966
+ ,127.65
+ ,52
+ ,537
+ ,1.319
+ ,0.89756
+ ,130.25
+ ,53
+ ,524
+ ,1.365
+ ,0.88444
+ ,131.85
+ ,54
+ ,536
+ ,14.016
+ ,0.8567
+ ,135.39
+ ,55
+ ,587
+ ,14.088
+ ,0.86092
+ ,133.09
+ ,56
+ ,597
+ ,14.268
+ ,0.86265
+ ,135.31
+ ,57
+ ,581
+ ,14.562
+ ,0.89135
+ ,133.14
+ ,58
+ ,564
+ ,14.816
+ ,0.91557
+ ,133.91
+ ,59
+ ,558
+ ,14.914
+ ,0.89892
+ ,132.97
+ ,60
+ ,575
+ ,14.614
+ ,0.89972
+ ,131.21
+ ,61
+ ,580
+ ,14.272
+ ,0.88305
+ ,130.34)
+ ,dim=c(5
+ ,61)
+ ,dimnames=list(c('Maand'
+ ,'Werkloosheidsgraad'
+ ,'Dollar/euro'
+ ,'Pond/euro'
+ ,'Yen/euro
')
+ ,1:61))
> y <- array(NA,dim=c(5,61),dimnames=list(c('Maand','Werkloosheidsgraad','Dollar/euro','Pond/euro','Yen/euro
'),1:61))
> for (i in 1:dim(x)[1])
+ {
+ for (j in 1:dim(x)[2])
+ {
+ y[i,j] <- as.numeric(x[i,j])
+ }
+ }
> par3 = 'No Linear Trend'
> par2 = 'Do not include Seasonal Dummies'
> par1 = '2'
> par3 <- 'No Linear Trend'
> par2 <- 'Do not include Seasonal Dummies'
> par1 <- '2'
> #'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, 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
Werkloosheidsgraad Maand Dollar/euro Pond/euro Yen/euro\r
1 591 1 13.119 0.69867 135.63
2 589 2 13.014 0.68968 136.55
3 584 3 13.201 0.69233 138.83
4 573 4 12.938 0.68293 138.84
5 567 5 12.694 0.68399 135.37
6 569 6 12.165 0.66895 132.22
7 621 7 12.037 0.68756 134.75
8 629 8 12.292 0.68527 135.98
9 628 9 12.256 0.67760 136.06
10 612 10 12.015 0.68137 138.05
11 595 11 11.786 0.67933 139.59
12 597 12 11.856 0.67922 140.58
13 593 13 12.103 0.68598 139.81
14 590 14 11.938 0.68297 140.77
15 580 15 1.202 0.68935 140.96
16 574 16 12.271 0.69463 143.59
17 573 17 1.277 0.68330 142.70
18 573 18 1.265 0.68666 145.11
19 620 19 12.684 0.68782 146.70
20 626 20 12.811 0.67669 148.53
21 620 21 12.727 0.67511 148.99
22 588 22 12.611 0.67254 149.65
23 566 23 12.881 0.67397 151.11
24 557 24 13.213 0.67286 154.82
25 561 25 12.999 0.66341 156.56
26 549 26 13.074 0.66800 157.60
27 532 27 13.242 0.68021 155.24
28 526 28 13.516 0.67934 160.68
29 511 29 13.511 0.68136 163.22
30 499 30 13.419 0.67562 164.55
31 555 31 13.716 0.67440 166.76
32 565 32 13.622 0.67766 159.05
33 542 33 13.896 0.68887 159.82
34 527 34 14.227 0.69614 164.95
35 510 35 14.684 0.70896 162.89
36 514 36 1.457 0.72064 163.55
37 517 37 14.718 0.74725 158.68
38 508 38 14.748 0.75094 157.97
39 493 39 15.527 0.77494 156.59
40 490 40 1.575 0.79487 161.56
41 469 41 15.557 0.79209 162.31
42 478 42 15.553 0.79152 166.26
43 528 43 1.577 0.79308 168.45
44 534 44 14.975 0.79279 163.63
45 518 45 14.369 0.79924 153.20
46 506 46 13.322 0.78668 133.52
47 502 47 12.732 0.83063 123.28
48 516 48 13.449 0.90448 122.51
49 528 49 13.239 0.91819 119.73
50 533 50 12.785 0.88691 118.30
51 536 51 1.305 0.91966 127.65
52 537 52 1.319 0.89756 130.25
53 524 53 1.365 0.88444 131.85
54 536 54 14.016 0.85670 135.39
55 587 55 14.088 0.86092 133.09
56 597 56 14.268 0.86265 135.31
57 581 57 14.562 0.89135 133.14
58 564 58 14.816 0.91557 133.91
59 558 59 14.914 0.89892 132.97
60 575 60 14.614 0.89972 131.21
61 580 61 14.272 0.88305 130.34
> k <- length(x[1,])
> df <- as.data.frame(x)
> (mylm <- lm(df))
Call:
lm(formula = df)
Coefficients:
(Intercept) Maand `Dollar/euro` `Pond/euro` `Yen/euro\\r`
1.092e+03 -4.732e-03 3.120e-01 -3.102e+02 -2.120e+00
> (mysum <- summary(mylm))
Call:
lm(formula = df)
Residuals:
Min 1Q Median 3Q Max
-75.195 -19.176 -2.138 21.132 54.799
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 1.092e+03 1.587e+02 6.883 5.40e-09 ***
Maand -4.732e-03 6.564e-01 -0.007 0.9943
`Dollar/euro` 3.120e-01 8.838e-01 0.353 0.7254
`Pond/euro` -3.102e+02 1.504e+02 -2.063 0.0438 *
`Yen/euro\\r` -2.120e+00 4.844e-01 -4.376 5.33e-05 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 29.59 on 56 degrees of freedom
Multiple R-squared: 0.5002, Adjusted R-squared: 0.4645
F-statistic: 14.01 on 4 and 56 DF, p-value: 5.541e-08
> 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.12749007 0.25498013 0.87250993
[2,] 0.10339569 0.20679138 0.89660431
[3,] 0.07828638 0.15657277 0.92171362
[4,] 0.05380350 0.10760699 0.94619650
[5,] 0.03183820 0.06367640 0.96816180
[6,] 0.07293021 0.14586042 0.92706979
[7,] 0.05181208 0.10362415 0.94818792
[8,] 0.03651113 0.07302227 0.96348887
[9,] 0.05290913 0.10581826 0.94709087
[10,] 0.02976309 0.05952619 0.97023691
[11,] 0.01629851 0.03259701 0.98370149
[12,] 0.02141207 0.04282413 0.97858793
[13,] 0.03928199 0.07856399 0.96071801
[14,] 0.06146424 0.12292849 0.93853576
[15,] 0.08966695 0.17933391 0.91033305
[16,] 0.15963784 0.31927568 0.84036216
[17,] 0.19940345 0.39880690 0.80059655
[18,] 0.20888740 0.41777480 0.79111260
[19,] 0.22511044 0.45022088 0.77488956
[20,] 0.41822585 0.83645170 0.58177415
[21,] 0.42242415 0.84484831 0.57757585
[22,] 0.38925114 0.77850228 0.61074886
[23,] 0.35982710 0.71965421 0.64017290
[24,] 0.51097996 0.97804008 0.48902004
[25,] 0.64572604 0.70854793 0.35427396
[26,] 0.73810828 0.52378344 0.26189172
[27,] 0.76089339 0.47821322 0.23910661
[28,] 0.76186296 0.47627408 0.23813704
[29,] 0.74316498 0.51367004 0.25683502
[30,] 0.79141330 0.41717340 0.20858670
[31,] 0.80363323 0.39273354 0.19636677
[32,] 0.76033308 0.47933383 0.23966692
[33,] 0.70859289 0.58281423 0.29140711
[34,] 0.73850838 0.52298325 0.26149162
[35,] 0.87119695 0.25760611 0.12880305
[36,] 0.90485075 0.19029851 0.09514925
[37,] 0.89744221 0.20511557 0.10255779
[38,] 0.84497733 0.31004533 0.15502267
[39,] 0.83636880 0.32726239 0.16363120
[40,] 0.84616635 0.30766730 0.15383365
[41,] 0.79537834 0.40924331 0.20462166
[42,] 0.70980535 0.58038929 0.29019465
[43,] 0.69076199 0.61847602 0.30923801
[44,] 0.57403722 0.85192555 0.42596278
[45,] 0.44872961 0.89745922 0.55127039
[46,] 0.29637175 0.59274350 0.70362825
> postscript(file="/var/wessaorg/rcomp/tmp/1ffiv1353881101.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/2yzv81353881101.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/32a9r1353881101.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/4dh7b1353881101.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/5oz2h1353881101.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 = 61
Frequency = 1
1 2 3 4 5 6
-1.28295435 -4.08359735 -3.48240122 -17.28983357 -30.23556004 -39.40749225
7 8 9 10 11 12
23.77191768 33.59407053 30.40075324 19.86814382 5.57595004 9.62322530
13 14 15 16 17 18
6.01534647 4.17291399 -0.09117955 -2.32761411 -5.29325584 0.86579670
19 20 21 22 23 24
48.03785460 54.43002786 48.94599055 17.58882728 -0.95240829 -2.53145377
25 26 27 28 29 30
2.29739350 -6.09319715 -24.35642303 -19.17586076 -28.15903389 -39.08666682
31 32 33 34 35 36
21.13154614 15.83383705 -2.13796137 -4.10766087 -21.63595467 -8.48279035
37 38 39 40 41 42
-11.68534294 -21.05049852 -31.77037648 -13.69637179 -38.32651854 -21.12452216
43 44 45 46 47 48
38.36674151 29.88436979 -6.02978721 -63.30951418 -75.19526172 -40.14175758
49 50 51 52 53 54
-29.71209242 -37.29839731 -0.73525881 -1.07801816 -14.76530317 -7.80758622
55 56 57 58 59 60
39.60821625 54.79907577 43.01366835 35.08317334 21.90080177 35.51658964
61
33.61367533
> postscript(file="/var/wessaorg/rcomp/tmp/6vmwb1353881101.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 = 61
Frequency = 1
lag(myerror, k = 1) myerror
0 -1.28295435 NA
1 -4.08359735 -1.28295435
2 -3.48240122 -4.08359735
3 -17.28983357 -3.48240122
4 -30.23556004 -17.28983357
5 -39.40749225 -30.23556004
6 23.77191768 -39.40749225
7 33.59407053 23.77191768
8 30.40075324 33.59407053
9 19.86814382 30.40075324
10 5.57595004 19.86814382
11 9.62322530 5.57595004
12 6.01534647 9.62322530
13 4.17291399 6.01534647
14 -0.09117955 4.17291399
15 -2.32761411 -0.09117955
16 -5.29325584 -2.32761411
17 0.86579670 -5.29325584
18 48.03785460 0.86579670
19 54.43002786 48.03785460
20 48.94599055 54.43002786
21 17.58882728 48.94599055
22 -0.95240829 17.58882728
23 -2.53145377 -0.95240829
24 2.29739350 -2.53145377
25 -6.09319715 2.29739350
26 -24.35642303 -6.09319715
27 -19.17586076 -24.35642303
28 -28.15903389 -19.17586076
29 -39.08666682 -28.15903389
30 21.13154614 -39.08666682
31 15.83383705 21.13154614
32 -2.13796137 15.83383705
33 -4.10766087 -2.13796137
34 -21.63595467 -4.10766087
35 -8.48279035 -21.63595467
36 -11.68534294 -8.48279035
37 -21.05049852 -11.68534294
38 -31.77037648 -21.05049852
39 -13.69637179 -31.77037648
40 -38.32651854 -13.69637179
41 -21.12452216 -38.32651854
42 38.36674151 -21.12452216
43 29.88436979 38.36674151
44 -6.02978721 29.88436979
45 -63.30951418 -6.02978721
46 -75.19526172 -63.30951418
47 -40.14175758 -75.19526172
48 -29.71209242 -40.14175758
49 -37.29839731 -29.71209242
50 -0.73525881 -37.29839731
51 -1.07801816 -0.73525881
52 -14.76530317 -1.07801816
53 -7.80758622 -14.76530317
54 39.60821625 -7.80758622
55 54.79907577 39.60821625
56 43.01366835 54.79907577
57 35.08317334 43.01366835
58 21.90080177 35.08317334
59 35.51658964 21.90080177
60 33.61367533 35.51658964
61 NA 33.61367533
> dum1 <- dum[2:length(myerror),]
> dum1
lag(myerror, k = 1) myerror
[1,] -4.08359735 -1.28295435
[2,] -3.48240122 -4.08359735
[3,] -17.28983357 -3.48240122
[4,] -30.23556004 -17.28983357
[5,] -39.40749225 -30.23556004
[6,] 23.77191768 -39.40749225
[7,] 33.59407053 23.77191768
[8,] 30.40075324 33.59407053
[9,] 19.86814382 30.40075324
[10,] 5.57595004 19.86814382
[11,] 9.62322530 5.57595004
[12,] 6.01534647 9.62322530
[13,] 4.17291399 6.01534647
[14,] -0.09117955 4.17291399
[15,] -2.32761411 -0.09117955
[16,] -5.29325584 -2.32761411
[17,] 0.86579670 -5.29325584
[18,] 48.03785460 0.86579670
[19,] 54.43002786 48.03785460
[20,] 48.94599055 54.43002786
[21,] 17.58882728 48.94599055
[22,] -0.95240829 17.58882728
[23,] -2.53145377 -0.95240829
[24,] 2.29739350 -2.53145377
[25,] -6.09319715 2.29739350
[26,] -24.35642303 -6.09319715
[27,] -19.17586076 -24.35642303
[28,] -28.15903389 -19.17586076
[29,] -39.08666682 -28.15903389
[30,] 21.13154614 -39.08666682
[31,] 15.83383705 21.13154614
[32,] -2.13796137 15.83383705
[33,] -4.10766087 -2.13796137
[34,] -21.63595467 -4.10766087
[35,] -8.48279035 -21.63595467
[36,] -11.68534294 -8.48279035
[37,] -21.05049852 -11.68534294
[38,] -31.77037648 -21.05049852
[39,] -13.69637179 -31.77037648
[40,] -38.32651854 -13.69637179
[41,] -21.12452216 -38.32651854
[42,] 38.36674151 -21.12452216
[43,] 29.88436979 38.36674151
[44,] -6.02978721 29.88436979
[45,] -63.30951418 -6.02978721
[46,] -75.19526172 -63.30951418
[47,] -40.14175758 -75.19526172
[48,] -29.71209242 -40.14175758
[49,] -37.29839731 -29.71209242
[50,] -0.73525881 -37.29839731
[51,] -1.07801816 -0.73525881
[52,] -14.76530317 -1.07801816
[53,] -7.80758622 -14.76530317
[54,] 39.60821625 -7.80758622
[55,] 54.79907577 39.60821625
[56,] 43.01366835 54.79907577
[57,] 35.08317334 43.01366835
[58,] 21.90080177 35.08317334
[59,] 35.51658964 21.90080177
[60,] 33.61367533 35.51658964
> z <- as.data.frame(dum1)
> z
lag(myerror, k = 1) myerror
1 -4.08359735 -1.28295435
2 -3.48240122 -4.08359735
3 -17.28983357 -3.48240122
4 -30.23556004 -17.28983357
5 -39.40749225 -30.23556004
6 23.77191768 -39.40749225
7 33.59407053 23.77191768
8 30.40075324 33.59407053
9 19.86814382 30.40075324
10 5.57595004 19.86814382
11 9.62322530 5.57595004
12 6.01534647 9.62322530
13 4.17291399 6.01534647
14 -0.09117955 4.17291399
15 -2.32761411 -0.09117955
16 -5.29325584 -2.32761411
17 0.86579670 -5.29325584
18 48.03785460 0.86579670
19 54.43002786 48.03785460
20 48.94599055 54.43002786
21 17.58882728 48.94599055
22 -0.95240829 17.58882728
23 -2.53145377 -0.95240829
24 2.29739350 -2.53145377
25 -6.09319715 2.29739350
26 -24.35642303 -6.09319715
27 -19.17586076 -24.35642303
28 -28.15903389 -19.17586076
29 -39.08666682 -28.15903389
30 21.13154614 -39.08666682
31 15.83383705 21.13154614
32 -2.13796137 15.83383705
33 -4.10766087 -2.13796137
34 -21.63595467 -4.10766087
35 -8.48279035 -21.63595467
36 -11.68534294 -8.48279035
37 -21.05049852 -11.68534294
38 -31.77037648 -21.05049852
39 -13.69637179 -31.77037648
40 -38.32651854 -13.69637179
41 -21.12452216 -38.32651854
42 38.36674151 -21.12452216
43 29.88436979 38.36674151
44 -6.02978721 29.88436979
45 -63.30951418 -6.02978721
46 -75.19526172 -63.30951418
47 -40.14175758 -75.19526172
48 -29.71209242 -40.14175758
49 -37.29839731 -29.71209242
50 -0.73525881 -37.29839731
51 -1.07801816 -0.73525881
52 -14.76530317 -1.07801816
53 -7.80758622 -14.76530317
54 39.60821625 -7.80758622
55 54.79907577 39.60821625
56 43.01366835 54.79907577
57 35.08317334 43.01366835
58 21.90080177 35.08317334
59 35.51658964 21.90080177
60 33.61367533 35.51658964
> 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/7zga51353881101.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/82deo1353881101.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/9cxcc1353881101.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/106iuo1353881101.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/115w071353881101.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/12y08a1353881101.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/13m2yc1353881101.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/14bq0m1353881101.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/15z0ba1353881101.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/16fftw1353881101.tab")
+ }
>
> try(system("convert tmp/1ffiv1353881101.ps tmp/1ffiv1353881101.png",intern=TRUE))
character(0)
> try(system("convert tmp/2yzv81353881101.ps tmp/2yzv81353881101.png",intern=TRUE))
character(0)
> try(system("convert tmp/32a9r1353881101.ps tmp/32a9r1353881101.png",intern=TRUE))
character(0)
> try(system("convert tmp/4dh7b1353881101.ps tmp/4dh7b1353881101.png",intern=TRUE))
character(0)
> try(system("convert tmp/5oz2h1353881101.ps tmp/5oz2h1353881101.png",intern=TRUE))
character(0)
> try(system("convert tmp/6vmwb1353881101.ps tmp/6vmwb1353881101.png",intern=TRUE))
character(0)
> try(system("convert tmp/7zga51353881101.ps tmp/7zga51353881101.png",intern=TRUE))
character(0)
> try(system("convert tmp/82deo1353881101.ps tmp/82deo1353881101.png",intern=TRUE))
character(0)
> try(system("convert tmp/9cxcc1353881101.ps tmp/9cxcc1353881101.png",intern=TRUE))
character(0)
> try(system("convert tmp/106iuo1353881101.ps tmp/106iuo1353881101.png",intern=TRUE))
character(0)
>
>
> proc.time()
user system elapsed
7.312 1.176 8.482