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(593530
+ ,3922
+ ,18004
+ ,707169
+ ,610763
+ ,3759
+ ,17537
+ ,703434
+ ,612613
+ ,4138
+ ,20366
+ ,701017
+ ,611324
+ ,4634
+ ,22782
+ ,696968
+ ,594167
+ ,3996
+ ,19169
+ ,688558
+ ,595454
+ ,4308
+ ,13807
+ ,679237
+ ,590865
+ ,4143
+ ,29743
+ ,677362
+ ,589379
+ ,4429
+ ,25591
+ ,676693
+ ,584428
+ ,5219
+ ,29096
+ ,670009
+ ,573100
+ ,4929
+ ,26482
+ ,667209
+ ,567456
+ ,5761
+ ,22405
+ ,662976
+ ,569028
+ ,5592
+ ,27044
+ ,660194
+ ,620735
+ ,4163
+ ,17970
+ ,652270
+ ,628884
+ ,4962
+ ,18730
+ ,648024
+ ,628232
+ ,5208
+ ,19684
+ ,629295
+ ,612117
+ ,4755
+ ,19785
+ ,624961
+ ,595404
+ ,4491
+ ,18479
+ ,617306
+ ,597141
+ ,5732
+ ,10698
+ ,607691
+ ,593408
+ ,5731
+ ,31956
+ ,596219
+ ,590072
+ ,5040
+ ,29506
+ ,591130
+ ,579799
+ ,6102
+ ,34506
+ ,584528
+ ,574205
+ ,4904
+ ,27165
+ ,576798
+ ,572775
+ ,5369
+ ,26736
+ ,575683
+ ,572942
+ ,5578
+ ,23691
+ ,574369
+ ,619567
+ ,4619
+ ,18157
+ ,566815
+ ,625809
+ ,4731
+ ,17328
+ ,573074
+ ,619916
+ ,5011
+ ,18205
+ ,567739
+ ,587625
+ ,5299
+ ,20995
+ ,571942
+ ,565742
+ ,4146
+ ,17382
+ ,570274
+ ,557274
+ ,4625
+ ,9367
+ ,568800
+ ,560576
+ ,4736
+ ,31124
+ ,558115
+ ,548854
+ ,4219
+ ,26551
+ ,550591
+ ,531673
+ ,5116
+ ,30651
+ ,548872
+ ,525919
+ ,4205
+ ,25859
+ ,547009
+ ,511038
+ ,4121
+ ,25100
+ ,545946
+ ,498662
+ ,5103
+ ,25778
+ ,539702
+ ,555362
+ ,4300
+ ,20418
+ ,542427
+ ,564591
+ ,4578
+ ,18688
+ ,542968
+ ,541657
+ ,3809
+ ,20424
+ ,536640
+ ,527070
+ ,5526
+ ,24776
+ ,533653
+ ,509846
+ ,4248
+ ,19814
+ ,540996
+ ,514258
+ ,3830
+ ,12738
+ ,538316
+ ,516922
+ ,4428
+ ,31566
+ ,532646
+ ,507561
+ ,4834
+ ,30111
+ ,533390
+ ,492622
+ ,4406
+ ,30019
+ ,528715
+ ,490243
+ ,4565
+ ,31934
+ ,530664
+ ,469357
+ ,4104
+ ,25826
+ ,528564
+ ,477580
+ ,4798
+ ,26835
+ ,519107
+ ,528379
+ ,3935
+ ,20205
+ ,518703
+ ,533590
+ ,3792
+ ,17789
+ ,519059
+ ,517945
+ ,4387
+ ,20520
+ ,518498
+ ,506174
+ ,4006
+ ,22518
+ ,524575
+ ,501866
+ ,4078
+ ,15572
+ ,536046
+ ,516141
+ ,4724
+ ,11509
+ ,552006
+ ,528222
+ ,3157
+ ,25447
+ ,560687
+ ,532638
+ ,3558
+ ,24090
+ ,578884
+ ,536322
+ ,3899
+ ,27786
+ ,591491
+ ,536535
+ ,4118
+ ,26195
+ ,599228
+ ,523597
+ ,3790
+ ,20516
+ ,633019
+ ,536214
+ ,4278
+ ,22759
+ ,649918
+ ,586570
+ ,4035
+ ,19028
+ ,655509)
+ ,dim=c(4
+ ,61)
+ ,dimnames=list(c('Werkzoekend'
+ ,'Bouw'
+ ,'Auto'
+ ,'Krediet')
+ ,1:61))
> y <- array(NA,dim=c(4,61),dimnames=list(c('Werkzoekend','Bouw','Auto','Krediet'),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 = '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
Werkzoekend Bouw Auto Krediet M1 M2 M3 M4 M5 M6 M7 M8 M9 M10 M11
1 593530 3922 18004 707169 1 0 0 0 0 0 0 0 0 0 0
2 610763 3759 17537 703434 0 1 0 0 0 0 0 0 0 0 0
3 612613 4138 20366 701017 0 0 1 0 0 0 0 0 0 0 0
4 611324 4634 22782 696968 0 0 0 1 0 0 0 0 0 0 0
5 594167 3996 19169 688558 0 0 0 0 1 0 0 0 0 0 0
6 595454 4308 13807 679237 0 0 0 0 0 1 0 0 0 0 0
7 590865 4143 29743 677362 0 0 0 0 0 0 1 0 0 0 0
8 589379 4429 25591 676693 0 0 0 0 0 0 0 1 0 0 0
9 584428 5219 29096 670009 0 0 0 0 0 0 0 0 1 0 0
10 573100 4929 26482 667209 0 0 0 0 0 0 0 0 0 1 0
11 567456 5761 22405 662976 0 0 0 0 0 0 0 0 0 0 1
12 569028 5592 27044 660194 0 0 0 0 0 0 0 0 0 0 0
13 620735 4163 17970 652270 1 0 0 0 0 0 0 0 0 0 0
14 628884 4962 18730 648024 0 1 0 0 0 0 0 0 0 0 0
15 628232 5208 19684 629295 0 0 1 0 0 0 0 0 0 0 0
16 612117 4755 19785 624961 0 0 0 1 0 0 0 0 0 0 0
17 595404 4491 18479 617306 0 0 0 0 1 0 0 0 0 0 0
18 597141 5732 10698 607691 0 0 0 0 0 1 0 0 0 0 0
19 593408 5731 31956 596219 0 0 0 0 0 0 1 0 0 0 0
20 590072 5040 29506 591130 0 0 0 0 0 0 0 1 0 0 0
21 579799 6102 34506 584528 0 0 0 0 0 0 0 0 1 0 0
22 574205 4904 27165 576798 0 0 0 0 0 0 0 0 0 1 0
23 572775 5369 26736 575683 0 0 0 0 0 0 0 0 0 0 1
24 572942 5578 23691 574369 0 0 0 0 0 0 0 0 0 0 0
25 619567 4619 18157 566815 1 0 0 0 0 0 0 0 0 0 0
26 625809 4731 17328 573074 0 1 0 0 0 0 0 0 0 0 0
27 619916 5011 18205 567739 0 0 1 0 0 0 0 0 0 0 0
28 587625 5299 20995 571942 0 0 0 1 0 0 0 0 0 0 0
29 565742 4146 17382 570274 0 0 0 0 1 0 0 0 0 0 0
30 557274 4625 9367 568800 0 0 0 0 0 1 0 0 0 0 0
31 560576 4736 31124 558115 0 0 0 0 0 0 1 0 0 0 0
32 548854 4219 26551 550591 0 0 0 0 0 0 0 1 0 0 0
33 531673 5116 30651 548872 0 0 0 0 0 0 0 0 1 0 0
34 525919 4205 25859 547009 0 0 0 0 0 0 0 0 0 1 0
35 511038 4121 25100 545946 0 0 0 0 0 0 0 0 0 0 1
36 498662 5103 25778 539702 0 0 0 0 0 0 0 0 0 0 0
37 555362 4300 20418 542427 1 0 0 0 0 0 0 0 0 0 0
38 564591 4578 18688 542968 0 1 0 0 0 0 0 0 0 0 0
39 541657 3809 20424 536640 0 0 1 0 0 0 0 0 0 0 0
40 527070 5526 24776 533653 0 0 0 1 0 0 0 0 0 0 0
41 509846 4248 19814 540996 0 0 0 0 1 0 0 0 0 0 0
42 514258 3830 12738 538316 0 0 0 0 0 1 0 0 0 0 0
43 516922 4428 31566 532646 0 0 0 0 0 0 1 0 0 0 0
44 507561 4834 30111 533390 0 0 0 0 0 0 0 1 0 0 0
45 492622 4406 30019 528715 0 0 0 0 0 0 0 0 1 0 0
46 490243 4565 31934 530664 0 0 0 0 0 0 0 0 0 1 0
47 469357 4104 25826 528564 0 0 0 0 0 0 0 0 0 0 1
48 477580 4798 26835 519107 0 0 0 0 0 0 0 0 0 0 0
49 528379 3935 20205 518703 1 0 0 0 0 0 0 0 0 0 0
50 533590 3792 17789 519059 0 1 0 0 0 0 0 0 0 0 0
51 517945 4387 20520 518498 0 0 1 0 0 0 0 0 0 0 0
52 506174 4006 22518 524575 0 0 0 1 0 0 0 0 0 0 0
53 501866 4078 15572 536046 0 0 0 0 1 0 0 0 0 0 0
54 516141 4724 11509 552006 0 0 0 0 0 1 0 0 0 0 0
55 528222 3157 25447 560687 0 0 0 0 0 0 1 0 0 0 0
56 532638 3558 24090 578884 0 0 0 0 0 0 0 1 0 0 0
57 536322 3899 27786 591491 0 0 0 0 0 0 0 0 1 0 0
58 536535 4118 26195 599228 0 0 0 0 0 0 0 0 0 1 0
59 523597 3790 20516 633019 0 0 0 0 0 0 0 0 0 0 1
60 536214 4278 22759 649918 0 0 0 0 0 0 0 0 0 0 0
61 586570 4035 19028 655509 1 0 0 0 0 0 0 0 0 0 0
> k <- length(x[1,])
> df <- as.data.frame(x)
> (mylm <- lm(df))
Call:
lm(formula = df)
Coefficients:
(Intercept) Bouw Auto Krediet M1 M2
2.094e+05 3.460e+01 -4.215e+00 4.288e-01 5.023e+04 5.216e+04
M3 M4 M5 M6 M7 M8
4.900e+04 3.218e+04 2.196e+04 -1.765e+04 7.051e+04 5.473e+04
M9 M10 M11
4.185e+04 3.894e+04 8.305e+03
> (mysum <- summary(mylm))
Call:
lm(formula = df)
Residuals:
Min 1Q Median 3Q Max
-39386 -10036 -2714 12879 35180
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 2.094e+05 5.025e+04 4.168 0.000134 ***
Bouw 3.460e+01 4.942e+00 7.001 9.09e-09 ***
Auto -4.215e+00 1.574e+00 -2.679 0.010217 *
Krediet 4.288e-01 4.561e-02 9.401 2.79e-12 ***
M1 5.023e+04 1.483e+04 3.387 0.001458 **
M2 5.216e+04 1.615e+04 3.229 0.002295 **
M3 4.900e+04 1.463e+04 3.349 0.001628 **
M4 3.218e+04 1.319e+04 2.440 0.018607 *
M5 2.196e+04 1.609e+04 1.364 0.179060
M6 -1.765e+04 2.407e+04 -0.733 0.467223
M7 7.051e+04 1.542e+04 4.572 3.64e-05 ***
M8 5.473e+04 1.348e+04 4.059 0.000189 ***
M9 4.185e+04 1.498e+04 2.793 0.007577 **
M10 3.894e+04 1.346e+04 2.892 0.005821 **
M11 8.305e+03 1.260e+04 0.659 0.513201
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 19630 on 46 degrees of freedom
Multiple R-squared: 0.8362, Adjusted R-squared: 0.7863
F-statistic: 16.77 on 14 and 46 DF, p-value: 1.456e-13
> 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.081771161 0.163542323 0.918228839
[2,] 0.091745402 0.183490803 0.908254598
[3,] 0.063370245 0.126740490 0.936629755
[4,] 0.034317564 0.068635128 0.965682436
[5,] 0.014827278 0.029654556 0.985172722
[6,] 0.008074494 0.016148989 0.991925506
[7,] 0.003674257 0.007348513 0.996325743
[8,] 0.002493072 0.004986145 0.997506928
[9,] 0.001566405 0.003132810 0.998433595
[10,] 0.001604411 0.003208821 0.998395589
[11,] 0.025143479 0.050286958 0.974856521
[12,] 0.183566416 0.367132832 0.816433584
[13,] 0.325891201 0.651782402 0.674108799
[14,] 0.456554867 0.913109733 0.543445133
[15,] 0.639641961 0.720716078 0.360358039
[16,] 0.713145785 0.573708429 0.286854215
[17,] 0.729152300 0.541695399 0.270847700
[18,] 0.899997431 0.200005137 0.100002569
[19,] 0.943571207 0.112857586 0.056428793
[20,] 0.975427520 0.049144961 0.024572480
[21,] 0.984710737 0.030578526 0.015289263
[22,] 0.986980537 0.026038926 0.013019463
[23,] 0.995596106 0.008807789 0.004403894
[24,] 0.992304136 0.015391728 0.007695864
[25,] 0.980229974 0.039540053 0.019770026
[26,] 0.955617590 0.088764820 0.044382410
> postscript(file="/var/www/html/rcomp/tmp/1i4mj1260641941.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/2e0ek1260641941.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/39wlz1260641941.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/4vb8d1260641941.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/5a2ox1260641941.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 = 61
Frequency = 1
1 2 3 4 5 6
-29107.4270 -8530.1582 -3674.8615 6619.9946 10129.4462 21622.6748
7 8 9 10 11 12
2565.0114 -10252.4547 -12015.5954 -20217.4351 -39385.7617 -2913.4408
13 14 15 16 17 18
13153.7361 -3246.5838 2798.5378 21466.2099 21880.4142 -8390.2357
19 20 21 22 23 24
-5717.6117 22488.3356 12258.8264 23395.4244 35179.7974 24148.5577
25 26 27 28 29 30
33635.4730 27895.9537 21456.4166 5984.4765 19696.1134 1108.5581
31 32 33 34 35 36
8707.0356 14601.5365 -2714.4039 6561.4483 22476.6853 -10035.7297
37 38 39 40 41 42
454.9719 -9387.4917 7473.7196 -30070.0915 -16924.4329 12879.1611
43 44 45 46 47 48
-11507.2404 -25588.7854 -11221.3909 -8954.4532 -8103.2515 -7279.1182
49 50 51 52 53 54
-4625.3256 -6731.7199 -28053.8124 -4000.5895 -34781.5409 -27220.1583
55 56 57 58 59 60
5952.8051 -1248.6320 13692.5639 -784.9844 -10167.4696 -3920.2690
61
-13511.4284
> postscript(file="/var/www/html/rcomp/tmp/6yf871260641941.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 = 61
Frequency = 1
lag(myerror, k = 1) myerror
0 -29107.4270 NA
1 -8530.1582 -29107.4270
2 -3674.8615 -8530.1582
3 6619.9946 -3674.8615
4 10129.4462 6619.9946
5 21622.6748 10129.4462
6 2565.0114 21622.6748
7 -10252.4547 2565.0114
8 -12015.5954 -10252.4547
9 -20217.4351 -12015.5954
10 -39385.7617 -20217.4351
11 -2913.4408 -39385.7617
12 13153.7361 -2913.4408
13 -3246.5838 13153.7361
14 2798.5378 -3246.5838
15 21466.2099 2798.5378
16 21880.4142 21466.2099
17 -8390.2357 21880.4142
18 -5717.6117 -8390.2357
19 22488.3356 -5717.6117
20 12258.8264 22488.3356
21 23395.4244 12258.8264
22 35179.7974 23395.4244
23 24148.5577 35179.7974
24 33635.4730 24148.5577
25 27895.9537 33635.4730
26 21456.4166 27895.9537
27 5984.4765 21456.4166
28 19696.1134 5984.4765
29 1108.5581 19696.1134
30 8707.0356 1108.5581
31 14601.5365 8707.0356
32 -2714.4039 14601.5365
33 6561.4483 -2714.4039
34 22476.6853 6561.4483
35 -10035.7297 22476.6853
36 454.9719 -10035.7297
37 -9387.4917 454.9719
38 7473.7196 -9387.4917
39 -30070.0915 7473.7196
40 -16924.4329 -30070.0915
41 12879.1611 -16924.4329
42 -11507.2404 12879.1611
43 -25588.7854 -11507.2404
44 -11221.3909 -25588.7854
45 -8954.4532 -11221.3909
46 -8103.2515 -8954.4532
47 -7279.1182 -8103.2515
48 -4625.3256 -7279.1182
49 -6731.7199 -4625.3256
50 -28053.8124 -6731.7199
51 -4000.5895 -28053.8124
52 -34781.5409 -4000.5895
53 -27220.1583 -34781.5409
54 5952.8051 -27220.1583
55 -1248.6320 5952.8051
56 13692.5639 -1248.6320
57 -784.9844 13692.5639
58 -10167.4696 -784.9844
59 -3920.2690 -10167.4696
60 -13511.4284 -3920.2690
61 NA -13511.4284
> dum1 <- dum[2:length(myerror),]
> dum1
lag(myerror, k = 1) myerror
[1,] -8530.1582 -29107.4270
[2,] -3674.8615 -8530.1582
[3,] 6619.9946 -3674.8615
[4,] 10129.4462 6619.9946
[5,] 21622.6748 10129.4462
[6,] 2565.0114 21622.6748
[7,] -10252.4547 2565.0114
[8,] -12015.5954 -10252.4547
[9,] -20217.4351 -12015.5954
[10,] -39385.7617 -20217.4351
[11,] -2913.4408 -39385.7617
[12,] 13153.7361 -2913.4408
[13,] -3246.5838 13153.7361
[14,] 2798.5378 -3246.5838
[15,] 21466.2099 2798.5378
[16,] 21880.4142 21466.2099
[17,] -8390.2357 21880.4142
[18,] -5717.6117 -8390.2357
[19,] 22488.3356 -5717.6117
[20,] 12258.8264 22488.3356
[21,] 23395.4244 12258.8264
[22,] 35179.7974 23395.4244
[23,] 24148.5577 35179.7974
[24,] 33635.4730 24148.5577
[25,] 27895.9537 33635.4730
[26,] 21456.4166 27895.9537
[27,] 5984.4765 21456.4166
[28,] 19696.1134 5984.4765
[29,] 1108.5581 19696.1134
[30,] 8707.0356 1108.5581
[31,] 14601.5365 8707.0356
[32,] -2714.4039 14601.5365
[33,] 6561.4483 -2714.4039
[34,] 22476.6853 6561.4483
[35,] -10035.7297 22476.6853
[36,] 454.9719 -10035.7297
[37,] -9387.4917 454.9719
[38,] 7473.7196 -9387.4917
[39,] -30070.0915 7473.7196
[40,] -16924.4329 -30070.0915
[41,] 12879.1611 -16924.4329
[42,] -11507.2404 12879.1611
[43,] -25588.7854 -11507.2404
[44,] -11221.3909 -25588.7854
[45,] -8954.4532 -11221.3909
[46,] -8103.2515 -8954.4532
[47,] -7279.1182 -8103.2515
[48,] -4625.3256 -7279.1182
[49,] -6731.7199 -4625.3256
[50,] -28053.8124 -6731.7199
[51,] -4000.5895 -28053.8124
[52,] -34781.5409 -4000.5895
[53,] -27220.1583 -34781.5409
[54,] 5952.8051 -27220.1583
[55,] -1248.6320 5952.8051
[56,] 13692.5639 -1248.6320
[57,] -784.9844 13692.5639
[58,] -10167.4696 -784.9844
[59,] -3920.2690 -10167.4696
[60,] -13511.4284 -3920.2690
> z <- as.data.frame(dum1)
> z
lag(myerror, k = 1) myerror
1 -8530.1582 -29107.4270
2 -3674.8615 -8530.1582
3 6619.9946 -3674.8615
4 10129.4462 6619.9946
5 21622.6748 10129.4462
6 2565.0114 21622.6748
7 -10252.4547 2565.0114
8 -12015.5954 -10252.4547
9 -20217.4351 -12015.5954
10 -39385.7617 -20217.4351
11 -2913.4408 -39385.7617
12 13153.7361 -2913.4408
13 -3246.5838 13153.7361
14 2798.5378 -3246.5838
15 21466.2099 2798.5378
16 21880.4142 21466.2099
17 -8390.2357 21880.4142
18 -5717.6117 -8390.2357
19 22488.3356 -5717.6117
20 12258.8264 22488.3356
21 23395.4244 12258.8264
22 35179.7974 23395.4244
23 24148.5577 35179.7974
24 33635.4730 24148.5577
25 27895.9537 33635.4730
26 21456.4166 27895.9537
27 5984.4765 21456.4166
28 19696.1134 5984.4765
29 1108.5581 19696.1134
30 8707.0356 1108.5581
31 14601.5365 8707.0356
32 -2714.4039 14601.5365
33 6561.4483 -2714.4039
34 22476.6853 6561.4483
35 -10035.7297 22476.6853
36 454.9719 -10035.7297
37 -9387.4917 454.9719
38 7473.7196 -9387.4917
39 -30070.0915 7473.7196
40 -16924.4329 -30070.0915
41 12879.1611 -16924.4329
42 -11507.2404 12879.1611
43 -25588.7854 -11507.2404
44 -11221.3909 -25588.7854
45 -8954.4532 -11221.3909
46 -8103.2515 -8954.4532
47 -7279.1182 -8103.2515
48 -4625.3256 -7279.1182
49 -6731.7199 -4625.3256
50 -28053.8124 -6731.7199
51 -4000.5895 -28053.8124
52 -34781.5409 -4000.5895
53 -27220.1583 -34781.5409
54 5952.8051 -27220.1583
55 -1248.6320 5952.8051
56 13692.5639 -1248.6320
57 -784.9844 13692.5639
58 -10167.4696 -784.9844
59 -3920.2690 -10167.4696
60 -13511.4284 -3920.2690
> 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/7d6a41260641941.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/85chw1260641941.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/95vv31260641941.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/10vwjf1260641941.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/11xcos1260641941.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/126e7r1260641941.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/1373ew1260641941.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/146l1k1260641941.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/15x52l1260641941.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/161jd31260641941.tab")
+ }
>
> try(system("convert tmp/1i4mj1260641941.ps tmp/1i4mj1260641941.png",intern=TRUE))
character(0)
> try(system("convert tmp/2e0ek1260641941.ps tmp/2e0ek1260641941.png",intern=TRUE))
character(0)
> try(system("convert tmp/39wlz1260641941.ps tmp/39wlz1260641941.png",intern=TRUE))
character(0)
> try(system("convert tmp/4vb8d1260641941.ps tmp/4vb8d1260641941.png",intern=TRUE))
character(0)
> try(system("convert tmp/5a2ox1260641941.ps tmp/5a2ox1260641941.png",intern=TRUE))
character(0)
> try(system("convert tmp/6yf871260641941.ps tmp/6yf871260641941.png",intern=TRUE))
character(0)
> try(system("convert tmp/7d6a41260641941.ps tmp/7d6a41260641941.png",intern=TRUE))
character(0)
> try(system("convert tmp/85chw1260641941.ps tmp/85chw1260641941.png",intern=TRUE))
character(0)
> try(system("convert tmp/95vv31260641941.ps tmp/95vv31260641941.png",intern=TRUE))
character(0)
> try(system("convert tmp/10vwjf1260641941.ps tmp/10vwjf1260641941.png",intern=TRUE))
character(0)
>
>
> proc.time()
user system elapsed
2.388 1.548 3.011