R version 2.12.0 (2010-10-15)
Copyright (C) 2010 The R Foundation for Statistical Computing
ISBN 3-900051-07-0
Platform: i486-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(99.2
+ ,96.7
+ ,101.0
+ ,99.0
+ ,98.1
+ ,100.1
+ ,631
+ ,923
+ ,-12
+ ,-10.8
+ ,654
+ ,294
+ ,-13
+ ,-12.2
+ ,671
+ ,833
+ ,-16
+ ,-14.1
+ ,586
+ ,840
+ ,-10
+ ,-15.2
+ ,600
+ ,969
+ ,-4
+ ,-15.8
+ ,625
+ ,568
+ ,-9
+ ,-15.8
+ ,558
+ ,110
+ ,-8
+ ,-14.9
+ ,630
+ ,577
+ ,-9
+ ,-12.6
+ ,628
+ ,654
+ ,-3
+ ,-9.9
+ ,603
+ ,184
+ ,-13
+ ,-7.8
+ ,656
+ ,255
+ ,-3
+ ,-6
+ ,600
+ ,730
+ ,-1
+ ,-5
+ ,670
+ ,326
+ ,-2
+ ,-4.5
+ ,678
+ ,423
+ ,0
+ ,-3.9
+ ,641
+ ,502
+ ,0
+ ,-2.9
+ ,625
+ ,311
+ ,-3
+ ,-1.5
+ ,628
+ ,177
+ ,0
+ ,-0.5
+ ,589
+ ,767
+ ,5
+ ,0
+ ,582
+ ,471
+ ,3
+ ,0.5
+ ,636
+ ,248
+ ,4
+ ,0.9
+ ,599
+ ,885
+ ,3
+ ,0.8
+ ,621
+ ,694
+ ,1
+ ,0.1
+ ,637
+ ,406
+ ,-1
+ ,-1
+ ,595
+ ,994
+ ,0
+ ,-2
+ ,696
+ ,308
+ ,-2
+ ,-3
+ ,674
+ ,201
+ ,-1
+ ,-3.7
+ ,648
+ ,861
+ ,2
+ ,-4.7
+ ,649
+ ,605
+ ,0
+ ,-6.4
+ ,672
+ ,392
+ ,-6
+ ,-7.5
+ ,598
+ ,396
+ ,-7
+ ,-7.8
+ ,613
+ ,177
+ ,-6
+ ,-7.7
+ ,638
+ ,104
+ ,-4
+ ,-6.6
+ ,615
+ ,632
+ ,-9
+ ,-4.2
+ ,634
+ ,465
+ ,-2
+ ,-2
+ ,638
+ ,686
+ ,-3
+ ,-0.7
+ ,604
+ ,243
+ ,2
+ ,0.1
+ ,706
+ ,669
+ ,3
+ ,0.9
+ ,677
+ ,185
+ ,1
+ ,2.1
+ ,644
+ ,328
+ ,0
+ ,3.5
+ ,644
+ ,825
+ ,1
+ ,4.9
+ ,605
+ ,707
+ ,1
+ ,5.7
+ ,600
+ ,136
+ ,3
+ ,6.2
+ ,612
+ ,166
+ ,5
+ ,6.5
+ ,599
+ ,659
+ ,5
+ ,6.5
+ ,634
+ ,210
+ ,4
+ ,6.3
+ ,618
+ ,234
+ ,11
+ ,6.2
+ ,613
+ ,576
+ ,8
+ ,6.4
+ ,627
+ ,200
+ ,-1
+ ,6.3
+ ,668
+ ,973
+ ,4
+ ,5.8
+ ,651
+ ,479
+ ,4
+ ,5.1
+ ,619
+ ,661
+ ,4
+ ,5.1
+ ,644
+ ,260
+ ,6
+ ,5.8
+ ,579
+ ,936
+ ,6
+ ,6.7
+ ,601
+ ,752
+ ,6
+ ,7.1
+ ,595
+ ,376
+ ,6
+ ,6.7
+ ,588
+ ,902
+ ,4
+ ,5.5
+ ,634
+ ,341
+ ,1
+ ,4.2
+ ,594
+ ,305
+ ,6
+ ,3
+ ,606
+ ,200
+ ,0
+ ,2.2
+ ,610
+ ,926
+ ,2
+ ,2
+ ,633
+ ,685
+ ,-2
+ ,1.8
+ ,639
+ ,696
+ ,0
+ ,1.8
+ ,659
+ ,451
+ ,1
+ ,1.5
+ ,593
+ ,248
+ ,-3
+ ,0.4
+ ,606
+ ,677
+ ,-3
+ ,-0.9
+ ,599
+ ,434
+ ,-5
+ ,-1.7
+ ,569
+ ,578
+ ,-7
+ ,-2.6
+ ,629
+ ,873
+ ,-7
+ ,-4.4
+ ,613
+ ,438
+ ,-5
+ ,-8.3
+ ,604
+ ,172
+ ,-13
+ ,-14.4
+ ,658
+ ,328
+ ,-16
+ ,-21.3
+ ,612
+ ,633
+ ,-20
+ ,-26.5
+ ,707
+ ,372
+ ,-18
+ ,-29.2
+ ,739
+ ,770
+ ,-21
+ ,-30.8
+ ,777
+ ,535
+ ,-20
+ ,-30.9
+ ,685
+ ,030
+ ,-16
+ ,-29.5
+ ,730
+ ,234
+ ,-14
+ ,-27.1
+ ,714
+ ,154
+ ,-12
+ ,-24.4
+ ,630
+ ,872
+ ,-10
+ ,-21.9
+ ,719
+ ,492
+ ,-3
+ ,-19.3
+ ,677
+ ,023
+ ,-4
+ ,-17
+ ,679
+ ,272
+ ,-4
+ ,-13.8
+ ,718
+ ,317
+ ,-1
+ ,-9.9
+ ,645
+ ,672
+ ,-8
+ ,-7.9)
+ ,dim=c(3
+ ,86)
+ ,dimnames=list(c('Werkloosheid'
+ ,'Consumenten'
+ ,'Ondernemers')
+ ,1:86))
> y <- array(NA,dim=c(3,86),dimnames=list(c('Werkloosheid','Consumenten','Ondernemers'),1:86))
> 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
> 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
Werkloosheid Consumenten Ondernemers M1 M2 M3 M4 M5 M6 M7 M8 M9 M10 M11
1 99.2 96.7 101.0 1 0 0 0 0 0 0 0 0 0 0
2 99.0 98.1 100.1 0 1 0 0 0 0 0 0 0 0 0
3 631.0 923.0 -12.0 0 0 1 0 0 0 0 0 0 0 0
4 -10.8 654.0 294.0 0 0 0 1 0 0 0 0 0 0 0
5 -13.0 -12.2 671.0 0 0 0 0 1 0 0 0 0 0 0
6 833.0 -16.0 -14.1 0 0 0 0 0 1 0 0 0 0 0
7 586.0 840.0 -10.0 0 0 0 0 0 0 1 0 0 0 0
8 -15.2 600.0 969.0 0 0 0 0 0 0 0 1 0 0 0
9 -4.0 -15.8 625.0 0 0 0 0 0 0 0 0 1 0 0
10 568.0 -9.0 -15.8 0 0 0 0 0 0 0 0 0 1 0
11 558.0 110.0 -8.0 0 0 0 0 0 0 0 0 0 0 1
12 -14.9 630.0 577.0 0 0 0 0 0 0 0 0 0 0 0
13 -9.0 -12.6 628.0 1 0 0 0 0 0 0 0 0 0 0
14 654.0 -3.0 -9.9 0 1 0 0 0 0 0 0 0 0 0
15 603.0 184.0 -13.0 0 0 1 0 0 0 0 0 0 0 0
16 -7.8 656.0 255.0 0 0 0 1 0 0 0 0 0 0 0
17 -3.0 -6.0 600.0 0 0 0 0 1 0 0 0 0 0 0
18 730.0 -1.0 -5.0 0 0 0 0 0 1 0 0 0 0 0
19 670.0 326.0 -2.0 0 0 0 0 0 0 1 0 0 0 0
20 -4.5 678.0 423.0 0 0 0 0 0 0 0 1 0 0 0
21 0.0 -3.9 641.0 0 0 0 0 0 0 0 0 1 0 0
22 502.0 0.0 -2.9 0 0 0 0 0 0 0 0 0 1 0
23 625.0 311.0 -3.0 0 0 0 0 0 0 0 0 0 0 1
24 -1.5 628.0 177.0 0 0 0 0 0 0 0 0 0 0 0
25 0.0 -0.5 589.0 1 0 0 0 0 0 0 0 0 0 0
26 767.0 5.0 0.0 0 1 0 0 0 0 0 0 0 0 0
27 582.0 471.0 3.0 0 0 1 0 0 0 0 0 0 0 0
28 0.5 636.0 248.0 0 0 0 1 0 0 0 0 0 0 0
29 4.0 0.9 599.0 0 0 0 0 1 0 0 0 0 0 0
30 885.0 3.0 0.8 0 0 0 0 0 1 0 0 0 0 0
31 621.0 694.0 1.0 0 0 0 0 0 0 1 0 0 0 0
32 0.1 637.0 406.0 0 0 0 0 0 0 0 1 0 0 0
33 -1.0 -1.0 595.0 0 0 0 0 0 0 0 0 1 0 0
34 994.0 0.0 -2.0 0 0 0 0 0 0 0 0 0 1 0
35 696.0 308.0 -2.0 0 0 0 0 0 0 0 0 0 0 1
36 -3.0 674.0 201.0 0 0 0 0 0 0 0 0 0 0 0
37 -1.0 -3.7 648.0 1 0 0 0 0 0 0 0 0 0 0
38 861.0 2.0 -4.7 0 1 0 0 0 0 0 0 0 0 0
39 649.0 605.0 0.0 0 0 1 0 0 0 0 0 0 0 0
40 -6.4 672.0 392.0 0 0 0 1 0 0 0 0 0 0 0
41 -6.0 -7.5 598.0 0 0 0 0 1 0 0 0 0 0 0
42 396.0 -7.0 -7.8 0 0 0 0 0 1 0 0 0 0 0
43 613.0 177.0 -6.0 0 0 0 0 0 0 1 0 0 0 0
44 -7.7 638.0 104.0 0 0 0 0 0 0 0 1 0 0 0
45 -4.0 -6.6 615.0 0 0 0 0 0 0 0 0 1 0 0
46 632.0 -9.0 -4.2 0 0 0 0 0 0 0 0 0 1 0
47 634.0 465.0 -2.0 0 0 0 0 0 0 0 0 0 0 1
48 -2.0 638.0 686.0 0 0 0 0 0 0 0 0 0 0 0
49 -3.0 -0.7 604.0 1 0 0 0 0 0 0 0 0 0 0
50 243.0 2.0 0.1 0 1 0 0 0 0 0 0 0 0 0
51 706.0 669.0 3.0 0 0 1 0 0 0 0 0 0 0 0
52 0.9 677.0 185.0 0 0 0 1 0 0 0 0 0 0 0
53 1.0 2.1 644.0 0 0 0 0 1 0 0 0 0 0 0
54 328.0 0.0 3.5 0 0 0 0 0 1 0 0 0 0 0
55 644.0 825.0 1.0 0 0 0 0 0 0 1 0 0 0 0
56 4.9 605.0 707.0 0 0 0 0 0 0 0 1 0 0 0
57 1.0 5.7 600.0 0 0 0 0 0 0 0 0 1 0 0
58 136.0 3.0 6.2 0 0 0 0 0 0 0 0 0 1 0
59 612.0 166.0 5.0 0 0 0 0 0 0 0 0 0 0 1
60 6.5 599.0 659.0 0 0 0 0 0 0 0 0 0 0 0
61 5.0 6.5 634.0 1 0 0 0 0 0 0 0 0 0 0
62 210.0 4.0 6.3 0 1 0 0 0 0 0 0 0 0 0
63 618.0 234.0 11.0 0 0 1 0 0 0 0 0 0 0 0
64 6.2 613.0 576.0 0 0 0 1 0 0 0 0 0 0 0
65 8.0 6.4 627.0 0 0 0 0 1 0 0 0 0 0 0
66 200.0 -1.0 6.3 0 0 0 0 0 1 0 0 0 0 0
67 668.0 973.0 4.0 0 0 0 0 0 0 1 0 0 0 0
68 5.8 651.0 479.0 0 0 0 0 0 0 0 1 0 0 0
69 4.0 5.1 619.0 0 0 0 0 0 0 0 0 1 0 0
70 661.0 4.0 5.1 0 0 0 0 0 0 0 0 0 1 0
71 644.0 260.0 6.0 0 0 0 0 0 0 0 0 0 0 1
72 5.8 579.0 936.0 0 0 0 0 0 0 0 0 0 0 0
73 6.0 6.7 601.0 1 0 0 0 0 0 0 0 0 0 0
74 752.0 6.0 7.1 0 1 0 0 0 0 0 0 0 0 0
75 595.0 376.0 6.0 0 0 1 0 0 0 0 0 0 0 0
76 6.7 588.0 902.0 0 0 0 1 0 0 0 0 0 0 0
77 4.0 5.5 634.0 0 0 0 0 1 0 0 0 0 0 0
78 341.0 1.0 4.2 0 0 0 0 0 1 0 0 0 0 0
79 594.0 305.0 6.0 0 0 0 0 0 0 1 0 0 0 0
80 3.0 606.0 200.0 0 0 0 0 0 0 0 1 0 0 0
81 0.0 2.2 610.0 0 0 0 0 0 0 0 0 1 0 0
82 926.0 2.0 2.0 0 0 0 0 0 0 0 0 0 1 0
83 633.0 685.0 -2.0 0 0 0 0 0 0 0 0 0 0 1
84 1.8 639.0 696.0 0 0 0 0 0 0 0 0 0 0 0
85 0.0 1.8 659.0 1 0 0 0 0 0 0 0 0 0 0
86 451.0 1.0 1.5 0 1 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) Consumenten Ondernemers M1 M2 M3
19.04415 0.01721 -0.05496 23.57176 486.02370 598.71434
M4 M5 M6 M7 M8 M9
-9.23348 14.60339 511.34105 598.73032 -6.02530 14.22094
M10 M11
612.17261 604.09891
> (mysum <- summary(mylm))
Call:
lm(formula = df)
Residuals:
Min 1Q Median 3Q Max
-494.93 -15.50 -3.24 14.38 362.67
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 19.04415 124.03495 0.154 0.878403
Consumenten 0.01721 0.13252 0.130 0.897034
Ondernemers -0.05496 0.11683 -0.470 0.639449
M1 23.57176 111.92681 0.211 0.833794
M2 486.02370 132.58290 3.666 0.000468 ***
M3 598.71434 105.40179 5.680 2.66e-07 ***
M4 -9.23348 81.20212 -0.114 0.909784
M5 14.60339 114.63019 0.127 0.898982
M6 511.34105 136.46072 3.747 0.000358 ***
M7 598.73032 103.30545 5.796 1.67e-07 ***
M8 -6.02530 79.91547 -0.075 0.940109
M9 14.22094 114.68974 0.124 0.901665
M10 612.17261 136.30688 4.491 2.64e-05 ***
M11 604.09891 112.43404 5.373 9.09e-07 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 148.2 on 72 degrees of freedom
Multiple R-squared: 0.8273, Adjusted R-squared: 0.7961
F-statistic: 26.53 on 13 and 72 DF, p-value: < 2.2e-16
> if (n > n25) {
+ kp3 <- k + 3
+ nmkm3 <- n - k - 3
+ gqarr <- array(NA, dim=c(nmkm3-kp3+1,3))
+ numgqtests <- 0
+ numsignificant1 <- 0
+ numsignificant5 <- 0
+ numsignificant10 <- 0
+ for (mypoint in kp3:nmkm3) {
+ j <- 0
+ numgqtests <- numgqtests + 1
+ for (myalt in c('greater', 'two.sided', 'less')) {
+ j <- j + 1
+ gqarr[mypoint-kp3+1,j] <- gqtest(mylm, point=mypoint, alternative=myalt)$p.value
+ }
+ if (gqarr[mypoint-kp3+1,2] < 0.01) numsignificant1 <- numsignificant1 + 1
+ if (gqarr[mypoint-kp3+1,2] < 0.05) numsignificant5 <- numsignificant5 + 1
+ if (gqarr[mypoint-kp3+1,2] < 0.10) numsignificant10 <- numsignificant10 + 1
+ }
+ gqarr
+ }
[,1] [,2] [,3]
[1,] 0.8707594 0.25848127 0.12924063
[2,] 0.8116996 0.37660080 0.18830040
[3,] 0.7069800 0.58604003 0.29302002
[4,] 0.6334488 0.73310236 0.36655118
[5,] 0.5137987 0.97240254 0.48620127
[6,] 0.4203941 0.84078816 0.57960592
[7,] 0.3311790 0.66235800 0.66882100
[8,] 0.2483141 0.49662828 0.75168586
[9,] 0.1735147 0.34702934 0.82648533
[10,] 0.4372154 0.87443073 0.56278463
[11,] 0.3525812 0.70516231 0.64741885
[12,] 0.2714492 0.54289838 0.72855081
[13,] 0.2021011 0.40420221 0.79789890
[14,] 0.3514922 0.70298436 0.64850782
[15,] 0.2763259 0.55265171 0.72367415
[16,] 0.2147823 0.42956462 0.78521769
[17,] 0.1591962 0.31839240 0.84080380
[18,] 0.6103725 0.77925501 0.38962750
[19,] 0.5606067 0.87878653 0.43939326
[20,] 0.4838809 0.96776183 0.51611908
[21,] 0.4073683 0.81473662 0.59263169
[22,] 0.7794758 0.44104846 0.22052423
[23,] 0.7229677 0.55406465 0.27703233
[24,] 0.6575145 0.68497102 0.34248551
[25,] 0.5860118 0.82797633 0.41398816
[26,] 0.7553130 0.48937397 0.24468699
[27,] 0.6956861 0.60862778 0.30431389
[28,] 0.6345836 0.73083282 0.36541641
[29,] 0.5616977 0.87660458 0.43830229
[30,] 0.4955481 0.99109624 0.50445188
[31,] 0.4209089 0.84181786 0.57909107
[32,] 0.3514192 0.70283836 0.64858082
[33,] 0.2836868 0.56737368 0.71631316
[34,] 0.4269417 0.85388332 0.57305834
[35,] 0.3677509 0.73550173 0.63224914
[36,] 0.2968324 0.59366487 0.70316757
[37,] 0.2324573 0.46491459 0.76754270
[38,] 0.2894208 0.57884160 0.71057920
[39,] 0.2248101 0.44962015 0.77518992
[40,] 0.1689926 0.33798511 0.83100745
[41,] 0.1218235 0.24364702 0.87817649
[42,] 0.8873678 0.22526434 0.11263217
[43,] 0.8356643 0.32867141 0.16433571
[44,] 0.7695127 0.46097451 0.23048725
[45,] 0.6878108 0.62437843 0.31218921
[46,] 0.9535606 0.09287872 0.04643936
[47,] 0.9195719 0.16085623 0.08042811
[48,] 0.8654141 0.26917184 0.13458592
[49,] 0.7861675 0.42766492 0.21383246
[50,] 0.7880830 0.42383403 0.21191701
[51,] 0.6772894 0.64542120 0.32271060
[52,] 0.5287164 0.94256715 0.47128358
[53,] 0.3610004 0.72200076 0.63899962
> postscript(file="/var/www/rcomp/tmp/1uf0y1292964344.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/www/rcomp/tmp/2mozj1292964344.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/www/rcomp/tmp/3mozj1292964344.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/www/rcomp/tmp/4mozj1292964344.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/www/rcomp/tmp/5xfym1292964344.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 = 86
Frequency = 1
1 2 3 4 5 6
60.4710150 -402.2544859 -3.3027624 -15.7072418 -9.5583707 302.1152042
7 8 9 10 11 12
-46.7803950 14.7129702 -2.6421968 -63.9302657 -67.4758403 -13.0735629
13 14 15 16 17 18
-16.8831993 148.4396630 -18.6396082 -14.8851626 -3.5673431 199.3572059
19 20 21 22 23 24
46.5051858 -5.9384154 2.0323905 -129.3761507 -3.6602223 -21.6237684
25 26 27 28 29 30
-10.2349400 261.8461033 -43.6994644 -6.6256955 3.2589471 354.6071434
31 32 33 34 35 36
-8.6631725 -1.5671561 -1.5457502 362.6733147 67.4463689 -22.5963462
37 38 39 40 41 42
-7.9371361 355.6394136 20.8295244 -6.2307868 -6.6514514 -134.6934270
43 44 45 46 47 48
-8.1503854 -25.9827581 -3.3501435 0.7072884 2.7444147 5.6795683
49 50 51 52 53 54
-12.4070746 -262.0967708 76.8929755 -10.3938799 2.7115655 -202.1928307
55 56 57 58 59 60
12.0823307 20.3269912 0.6137513 -494.9276301 -13.7250944 13.3667921
61 62 63 64 65 66
-2.8821390 -294.7904290 -3.1810257 17.4975248 8.7032164 -330.0217284
67 68 69 70 71 72
33.7001502 7.9040994 4.6683470 29.9947022 16.7121366 28.2353431
73 74 75 76 77 78
-3.6993125 247.2191205 -28.8996392 36.3452419 5.1034362 -189.1715675
79 80 81 82 83 84
-28.6937138 -9.4557312 0.2236016 294.8587412 -2.0417632 10.0119740
85 86
-6.4272134 -54.0026148
> postscript(file="/var/www/rcomp/tmp/6xfym1292964344.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 = 86
Frequency = 1
lag(myerror, k = 1) myerror
0 60.4710150 NA
1 -402.2544859 60.4710150
2 -3.3027624 -402.2544859
3 -15.7072418 -3.3027624
4 -9.5583707 -15.7072418
5 302.1152042 -9.5583707
6 -46.7803950 302.1152042
7 14.7129702 -46.7803950
8 -2.6421968 14.7129702
9 -63.9302657 -2.6421968
10 -67.4758403 -63.9302657
11 -13.0735629 -67.4758403
12 -16.8831993 -13.0735629
13 148.4396630 -16.8831993
14 -18.6396082 148.4396630
15 -14.8851626 -18.6396082
16 -3.5673431 -14.8851626
17 199.3572059 -3.5673431
18 46.5051858 199.3572059
19 -5.9384154 46.5051858
20 2.0323905 -5.9384154
21 -129.3761507 2.0323905
22 -3.6602223 -129.3761507
23 -21.6237684 -3.6602223
24 -10.2349400 -21.6237684
25 261.8461033 -10.2349400
26 -43.6994644 261.8461033
27 -6.6256955 -43.6994644
28 3.2589471 -6.6256955
29 354.6071434 3.2589471
30 -8.6631725 354.6071434
31 -1.5671561 -8.6631725
32 -1.5457502 -1.5671561
33 362.6733147 -1.5457502
34 67.4463689 362.6733147
35 -22.5963462 67.4463689
36 -7.9371361 -22.5963462
37 355.6394136 -7.9371361
38 20.8295244 355.6394136
39 -6.2307868 20.8295244
40 -6.6514514 -6.2307868
41 -134.6934270 -6.6514514
42 -8.1503854 -134.6934270
43 -25.9827581 -8.1503854
44 -3.3501435 -25.9827581
45 0.7072884 -3.3501435
46 2.7444147 0.7072884
47 5.6795683 2.7444147
48 -12.4070746 5.6795683
49 -262.0967708 -12.4070746
50 76.8929755 -262.0967708
51 -10.3938799 76.8929755
52 2.7115655 -10.3938799
53 -202.1928307 2.7115655
54 12.0823307 -202.1928307
55 20.3269912 12.0823307
56 0.6137513 20.3269912
57 -494.9276301 0.6137513
58 -13.7250944 -494.9276301
59 13.3667921 -13.7250944
60 -2.8821390 13.3667921
61 -294.7904290 -2.8821390
62 -3.1810257 -294.7904290
63 17.4975248 -3.1810257
64 8.7032164 17.4975248
65 -330.0217284 8.7032164
66 33.7001502 -330.0217284
67 7.9040994 33.7001502
68 4.6683470 7.9040994
69 29.9947022 4.6683470
70 16.7121366 29.9947022
71 28.2353431 16.7121366
72 -3.6993125 28.2353431
73 247.2191205 -3.6993125
74 -28.8996392 247.2191205
75 36.3452419 -28.8996392
76 5.1034362 36.3452419
77 -189.1715675 5.1034362
78 -28.6937138 -189.1715675
79 -9.4557312 -28.6937138
80 0.2236016 -9.4557312
81 294.8587412 0.2236016
82 -2.0417632 294.8587412
83 10.0119740 -2.0417632
84 -6.4272134 10.0119740
85 -54.0026148 -6.4272134
86 NA -54.0026148
> dum1 <- dum[2:length(myerror),]
> dum1
lag(myerror, k = 1) myerror
[1,] -402.2544859 60.4710150
[2,] -3.3027624 -402.2544859
[3,] -15.7072418 -3.3027624
[4,] -9.5583707 -15.7072418
[5,] 302.1152042 -9.5583707
[6,] -46.7803950 302.1152042
[7,] 14.7129702 -46.7803950
[8,] -2.6421968 14.7129702
[9,] -63.9302657 -2.6421968
[10,] -67.4758403 -63.9302657
[11,] -13.0735629 -67.4758403
[12,] -16.8831993 -13.0735629
[13,] 148.4396630 -16.8831993
[14,] -18.6396082 148.4396630
[15,] -14.8851626 -18.6396082
[16,] -3.5673431 -14.8851626
[17,] 199.3572059 -3.5673431
[18,] 46.5051858 199.3572059
[19,] -5.9384154 46.5051858
[20,] 2.0323905 -5.9384154
[21,] -129.3761507 2.0323905
[22,] -3.6602223 -129.3761507
[23,] -21.6237684 -3.6602223
[24,] -10.2349400 -21.6237684
[25,] 261.8461033 -10.2349400
[26,] -43.6994644 261.8461033
[27,] -6.6256955 -43.6994644
[28,] 3.2589471 -6.6256955
[29,] 354.6071434 3.2589471
[30,] -8.6631725 354.6071434
[31,] -1.5671561 -8.6631725
[32,] -1.5457502 -1.5671561
[33,] 362.6733147 -1.5457502
[34,] 67.4463689 362.6733147
[35,] -22.5963462 67.4463689
[36,] -7.9371361 -22.5963462
[37,] 355.6394136 -7.9371361
[38,] 20.8295244 355.6394136
[39,] -6.2307868 20.8295244
[40,] -6.6514514 -6.2307868
[41,] -134.6934270 -6.6514514
[42,] -8.1503854 -134.6934270
[43,] -25.9827581 -8.1503854
[44,] -3.3501435 -25.9827581
[45,] 0.7072884 -3.3501435
[46,] 2.7444147 0.7072884
[47,] 5.6795683 2.7444147
[48,] -12.4070746 5.6795683
[49,] -262.0967708 -12.4070746
[50,] 76.8929755 -262.0967708
[51,] -10.3938799 76.8929755
[52,] 2.7115655 -10.3938799
[53,] -202.1928307 2.7115655
[54,] 12.0823307 -202.1928307
[55,] 20.3269912 12.0823307
[56,] 0.6137513 20.3269912
[57,] -494.9276301 0.6137513
[58,] -13.7250944 -494.9276301
[59,] 13.3667921 -13.7250944
[60,] -2.8821390 13.3667921
[61,] -294.7904290 -2.8821390
[62,] -3.1810257 -294.7904290
[63,] 17.4975248 -3.1810257
[64,] 8.7032164 17.4975248
[65,] -330.0217284 8.7032164
[66,] 33.7001502 -330.0217284
[67,] 7.9040994 33.7001502
[68,] 4.6683470 7.9040994
[69,] 29.9947022 4.6683470
[70,] 16.7121366 29.9947022
[71,] 28.2353431 16.7121366
[72,] -3.6993125 28.2353431
[73,] 247.2191205 -3.6993125
[74,] -28.8996392 247.2191205
[75,] 36.3452419 -28.8996392
[76,] 5.1034362 36.3452419
[77,] -189.1715675 5.1034362
[78,] -28.6937138 -189.1715675
[79,] -9.4557312 -28.6937138
[80,] 0.2236016 -9.4557312
[81,] 294.8587412 0.2236016
[82,] -2.0417632 294.8587412
[83,] 10.0119740 -2.0417632
[84,] -6.4272134 10.0119740
[85,] -54.0026148 -6.4272134
> z <- as.data.frame(dum1)
> z
lag(myerror, k = 1) myerror
1 -402.2544859 60.4710150
2 -3.3027624 -402.2544859
3 -15.7072418 -3.3027624
4 -9.5583707 -15.7072418
5 302.1152042 -9.5583707
6 -46.7803950 302.1152042
7 14.7129702 -46.7803950
8 -2.6421968 14.7129702
9 -63.9302657 -2.6421968
10 -67.4758403 -63.9302657
11 -13.0735629 -67.4758403
12 -16.8831993 -13.0735629
13 148.4396630 -16.8831993
14 -18.6396082 148.4396630
15 -14.8851626 -18.6396082
16 -3.5673431 -14.8851626
17 199.3572059 -3.5673431
18 46.5051858 199.3572059
19 -5.9384154 46.5051858
20 2.0323905 -5.9384154
21 -129.3761507 2.0323905
22 -3.6602223 -129.3761507
23 -21.6237684 -3.6602223
24 -10.2349400 -21.6237684
25 261.8461033 -10.2349400
26 -43.6994644 261.8461033
27 -6.6256955 -43.6994644
28 3.2589471 -6.6256955
29 354.6071434 3.2589471
30 -8.6631725 354.6071434
31 -1.5671561 -8.6631725
32 -1.5457502 -1.5671561
33 362.6733147 -1.5457502
34 67.4463689 362.6733147
35 -22.5963462 67.4463689
36 -7.9371361 -22.5963462
37 355.6394136 -7.9371361
38 20.8295244 355.6394136
39 -6.2307868 20.8295244
40 -6.6514514 -6.2307868
41 -134.6934270 -6.6514514
42 -8.1503854 -134.6934270
43 -25.9827581 -8.1503854
44 -3.3501435 -25.9827581
45 0.7072884 -3.3501435
46 2.7444147 0.7072884
47 5.6795683 2.7444147
48 -12.4070746 5.6795683
49 -262.0967708 -12.4070746
50 76.8929755 -262.0967708
51 -10.3938799 76.8929755
52 2.7115655 -10.3938799
53 -202.1928307 2.7115655
54 12.0823307 -202.1928307
55 20.3269912 12.0823307
56 0.6137513 20.3269912
57 -494.9276301 0.6137513
58 -13.7250944 -494.9276301
59 13.3667921 -13.7250944
60 -2.8821390 13.3667921
61 -294.7904290 -2.8821390
62 -3.1810257 -294.7904290
63 17.4975248 -3.1810257
64 8.7032164 17.4975248
65 -330.0217284 8.7032164
66 33.7001502 -330.0217284
67 7.9040994 33.7001502
68 4.6683470 7.9040994
69 29.9947022 4.6683470
70 16.7121366 29.9947022
71 28.2353431 16.7121366
72 -3.6993125 28.2353431
73 247.2191205 -3.6993125
74 -28.8996392 247.2191205
75 36.3452419 -28.8996392
76 5.1034362 36.3452419
77 -189.1715675 5.1034362
78 -28.6937138 -189.1715675
79 -9.4557312 -28.6937138
80 0.2236016 -9.4557312
81 294.8587412 0.2236016
82 -2.0417632 294.8587412
83 10.0119740 -2.0417632
84 -6.4272134 10.0119740
85 -54.0026148 -6.4272134
> 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/rcomp/tmp/78of61292964344.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/www/rcomp/tmp/88of61292964344.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/www/rcomp/tmp/91xxr1292964344.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/www/rcomp/tmp/101xxr1292964344.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/www/rcomp/createtable file can be downloaded at http://www.wessa.net/cretab
> load(file="/var/www/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/rcomp/tmp/11mgdf1292964344.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/rcomp/tmp/127hu31292964344.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/rcomp/tmp/13e09x1292964344.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/rcomp/tmp/1479qi1292964344.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/rcomp/tmp/15s97o1292964344.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/rcomp/tmp/1661mw1292964344.tab")
+ }
>
> try(system("convert tmp/1uf0y1292964344.ps tmp/1uf0y1292964344.png",intern=TRUE))
character(0)
> try(system("convert tmp/2mozj1292964344.ps tmp/2mozj1292964344.png",intern=TRUE))
character(0)
> try(system("convert tmp/3mozj1292964344.ps tmp/3mozj1292964344.png",intern=TRUE))
character(0)
> try(system("convert tmp/4mozj1292964344.ps tmp/4mozj1292964344.png",intern=TRUE))
character(0)
> try(system("convert tmp/5xfym1292964344.ps tmp/5xfym1292964344.png",intern=TRUE))
character(0)
> try(system("convert tmp/6xfym1292964344.ps tmp/6xfym1292964344.png",intern=TRUE))
character(0)
> try(system("convert tmp/78of61292964344.ps tmp/78of61292964344.png",intern=TRUE))
character(0)
> try(system("convert tmp/88of61292964344.ps tmp/88of61292964344.png",intern=TRUE))
character(0)
> try(system("convert tmp/91xxr1292964344.ps tmp/91xxr1292964344.png",intern=TRUE))
character(0)
> try(system("convert tmp/101xxr1292964344.ps tmp/101xxr1292964344.png",intern=TRUE))
character(0)
>
>
> proc.time()
user system elapsed
3.560 1.450 4.989