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(96.8,610763,114.1,612613,110.3,611324,103.9,594167,101.6,595454,94.6,590865,95.9,589379,104.7,584428,102.8,573100,98.1,567456,113.9,569028,80.9,620735,95.7,628884,113.2,628232,105.9,612117,108.8,595404,102.3,597141,99,593408,100.7,590072,115.5,579799,100.7,574205,109.9,572775,114.6,572942,85.4,619567,100.5,625809,114.8,619916,116.5,587625,112.9,565742,102,557274,106,560576,105.3,548854,118.8,531673,106.1,525919,109.3,511038,117.2,498662,92.5,555362,104.2,564591,112.5,541657,122.4,527070,113.3,509846,100,514258,110.7,516922,112.8,507561,109.8,492622,117.3,490243,109.1,469357,115.9,477580,96,528379,99.8,533590,116.8,517945,115.7,506174,99.4,501866,94.3,516141,91,528222,93.2,532638,103.1,536322,94.1,536535,91.8,523597,102.7,536214,82.6,586570,89.1,596594),dim=c(2,61),dimnames=list(c('Tot_ind_productie','Tot_nietwerkende_werkzoekenden'),1:61))
> y <- array(NA,dim=c(2,61),dimnames=list(c('Tot_ind_productie','Tot_nietwerkende_werkzoekenden'),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 = '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.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
Tot_nietwerkende_werkzoekenden Tot_ind_productie M1 M2 M3 M4 M5 M6 M7 M8 M9
1 610763 96.8 1 0 0 0 0 0 0 0 0
2 612613 114.1 0 1 0 0 0 0 0 0 0
3 611324 110.3 0 0 1 0 0 0 0 0 0
4 594167 103.9 0 0 0 1 0 0 0 0 0
5 595454 101.6 0 0 0 0 1 0 0 0 0
6 590865 94.6 0 0 0 0 0 1 0 0 0
7 589379 95.9 0 0 0 0 0 0 1 0 0
8 584428 104.7 0 0 0 0 0 0 0 1 0
9 573100 102.8 0 0 0 0 0 0 0 0 1
10 567456 98.1 0 0 0 0 0 0 0 0 0
11 569028 113.9 0 0 0 0 0 0 0 0 0
12 620735 80.9 0 0 0 0 0 0 0 0 0
13 628884 95.7 1 0 0 0 0 0 0 0 0
14 628232 113.2 0 1 0 0 0 0 0 0 0
15 612117 105.9 0 0 1 0 0 0 0 0 0
16 595404 108.8 0 0 0 1 0 0 0 0 0
17 597141 102.3 0 0 0 0 1 0 0 0 0
18 593408 99.0 0 0 0 0 0 1 0 0 0
19 590072 100.7 0 0 0 0 0 0 1 0 0
20 579799 115.5 0 0 0 0 0 0 0 1 0
21 574205 100.7 0 0 0 0 0 0 0 0 1
22 572775 109.9 0 0 0 0 0 0 0 0 0
23 572942 114.6 0 0 0 0 0 0 0 0 0
24 619567 85.4 0 0 0 0 0 0 0 0 0
25 625809 100.5 1 0 0 0 0 0 0 0 0
26 619916 114.8 0 1 0 0 0 0 0 0 0
27 587625 116.5 0 0 1 0 0 0 0 0 0
28 565742 112.9 0 0 0 1 0 0 0 0 0
29 557274 102.0 0 0 0 0 1 0 0 0 0
30 560576 106.0 0 0 0 0 0 1 0 0 0
31 548854 105.3 0 0 0 0 0 0 1 0 0
32 531673 118.8 0 0 0 0 0 0 0 1 0
33 525919 106.1 0 0 0 0 0 0 0 0 1
34 511038 109.3 0 0 0 0 0 0 0 0 0
35 498662 117.2 0 0 0 0 0 0 0 0 0
36 555362 92.5 0 0 0 0 0 0 0 0 0
37 564591 104.2 1 0 0 0 0 0 0 0 0
38 541657 112.5 0 1 0 0 0 0 0 0 0
39 527070 122.4 0 0 1 0 0 0 0 0 0
40 509846 113.3 0 0 0 1 0 0 0 0 0
41 514258 100.0 0 0 0 0 1 0 0 0 0
42 516922 110.7 0 0 0 0 0 1 0 0 0
43 507561 112.8 0 0 0 0 0 0 1 0 0
44 492622 109.8 0 0 0 0 0 0 0 1 0
45 490243 117.3 0 0 0 0 0 0 0 0 1
46 469357 109.1 0 0 0 0 0 0 0 0 0
47 477580 115.9 0 0 0 0 0 0 0 0 0
48 528379 96.0 0 0 0 0 0 0 0 0 0
49 533590 99.8 1 0 0 0 0 0 0 0 0
50 517945 116.8 0 1 0 0 0 0 0 0 0
51 506174 115.7 0 0 1 0 0 0 0 0 0
52 501866 99.4 0 0 0 1 0 0 0 0 0
53 516141 94.3 0 0 0 0 1 0 0 0 0
54 528222 91.0 0 0 0 0 0 1 0 0 0
55 532638 93.2 0 0 0 0 0 0 1 0 0
56 536322 103.1 0 0 0 0 0 0 0 1 0
57 536535 94.1 0 0 0 0 0 0 0 0 1
58 523597 91.8 0 0 0 0 0 0 0 0 0
59 536214 102.7 0 0 0 0 0 0 0 0 0
60 586570 82.6 0 0 0 0 0 0 0 0 0
61 596594 89.1 1 0 0 0 0 0 0 0 0
M10 M11
1 0 0
2 0 0
3 0 0
4 0 0
5 0 0
6 0 0
7 0 0
8 0 0
9 0 0
10 1 0
11 0 1
12 0 0
13 0 0
14 0 0
15 0 0
16 0 0
17 0 0
18 0 0
19 0 0
20 0 0
21 0 0
22 1 0
23 0 1
24 0 0
25 0 0
26 0 0
27 0 0
28 0 0
29 0 0
30 0 0
31 0 0
32 0 0
33 0 0
34 1 0
35 0 1
36 0 0
37 0 0
38 0 0
39 0 0
40 0 0
41 0 0
42 0 0
43 0 0
44 0 0
45 0 0
46 1 0
47 0 1
48 0 0
49 0 0
50 0 0
51 0 0
52 0 0
53 0 0
54 0 0
55 0 0
56 0 0
57 0 0
58 1 0
59 0 1
60 0 0
61 0 0
> k <- length(x[1,])
> df <- as.data.frame(x)
> (mylm <- lm(df))
Call:
lm(formula = df)
Coefficients:
(Intercept) Tot_ind_productie M1 M2
760951.9 -2044.2 32107.2 56735.4
M3 M4 M5 M6
41279.5 12535.0 -393.5 2001.3
M7 M8 M9 M10
401.9 9659.1 -7942.7 -20243.2
M11
645.2
> (mysum <- summary(mylm))
Call:
lm(formula = df)
Residuals:
Min 1Q Median 3Q Max
-68424 -24948 3917 31458 56727
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 760951.9 78419.3 9.704 6.72e-13 ***
Tot_ind_productie -2044.2 873.4 -2.341 0.0235 *
M1 32107.2 25524.5 1.258 0.2145
M2 56735.4 34234.0 1.657 0.1040
M3 41279.5 34162.4 1.208 0.2328
M4 12535.0 30573.3 0.410 0.6836
M5 -393.5 27284.2 -0.014 0.9886
M6 2001.3 27362.0 0.073 0.9420
M7 401.9 27852.2 0.014 0.9885
M8 9659.1 32001.8 0.302 0.7641
M9 -7942.7 28936.8 -0.274 0.7849
M10 -20243.2 28693.1 -0.706 0.4839
M11 645.2 33398.3 0.019 0.9847
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 39500 on 48 degrees of freedom
Multiple R-squared: 0.3094, Adjusted R-squared: 0.1368
F-statistic: 1.792 on 12 and 48 DF, p-value: 0.07667
> 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,] 1.941949e-02 3.883899e-02 0.98058051
[2,] 4.449830e-03 8.899659e-03 0.99555017
[3,] 9.865369e-04 1.973074e-03 0.99901346
[4,] 2.021813e-04 4.043625e-04 0.99979782
[5,] 4.336201e-05 8.672402e-05 0.99995664
[6,] 8.283140e-06 1.656628e-05 0.99999172
[7,] 3.205506e-06 6.411012e-06 0.99999679
[8,] 1.055813e-06 2.111626e-06 0.99999894
[9,] 2.714781e-07 5.429562e-07 0.99999973
[10,] 1.332747e-07 2.665493e-07 0.99999987
[11,] 1.273109e-07 2.546218e-07 0.99999987
[12,] 6.255319e-06 1.251064e-05 0.99999374
[13,] 1.701019e-04 3.402038e-04 0.99982990
[14,] 7.694674e-03 1.538935e-02 0.99230533
[15,] 1.814121e-02 3.628241e-02 0.98185879
[16,] 4.441661e-02 8.883322e-02 0.95558339
[17,] 1.206305e-01 2.412609e-01 0.87936955
[18,] 1.809586e-01 3.619173e-01 0.81904136
[19,] 3.561489e-01 7.122978e-01 0.64385108
[20,] 5.561217e-01 8.877567e-01 0.44387834
[21,] 5.197101e-01 9.605797e-01 0.48028987
[22,] 5.216900e-01 9.566200e-01 0.47831002
[23,] 7.246841e-01 5.506318e-01 0.27531591
[24,] 7.442088e-01 5.115824e-01 0.25579121
[25,] 8.298588e-01 3.402824e-01 0.17014122
[26,] 8.380352e-01 3.239295e-01 0.16196477
[27,] 8.718816e-01 2.562369e-01 0.12811843
[28,] 8.906189e-01 2.187622e-01 0.10938112
[29,] 8.879515e-01 2.240970e-01 0.11204850
[30,] 9.454170e-01 1.091660e-01 0.05458302
> postscript(file="/var/www/html/rcomp/tmp/1ymex1258623606.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/2rfxb1258623606.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/33nob1258623606.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/45wgy1258623606.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/52zbc1258623606.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 7
15585.429 28172.438 34571.269 33075.692 42589.400 21296.054 24066.969
8 9 10 11 12 13 14
27847.969 30237.677 27286.361 40268.800 25161.361 31457.776 41951.631
15 16 17 18 19 20 21
26369.654 44329.423 45707.362 32833.669 34572.277 45296.662 27049.792
22 23 24 25 26 27 28
56727.285 45613.762 33192.400 38195.083 36906.400 23546.500 23048.769
29 30 31 32 33 34 35
5227.092 14311.285 2757.739 3916.623 -10197.361 -6236.254 -23351.238
36 37 38 39 40 41 42
-16498.561 -15459.263 -46054.331 -24947.538 -32029.538 -41877.369 -19734.830
43 44 45 46 47 48 49
-23203.530 -53532.454 -22977.976 -48326.100 -47090.738 -36326.754 -55454.878
50 51 52 53 54 55 56
-60976.138 -59539.885 -68424.346 -51646.485 -48706.177 -38193.454 -23528.800
57 58 59 60 61
-24112.131 -29451.293 -15440.585 -5528.446 -14324.148
> postscript(file="/var/www/html/rcomp/tmp/63dl51258623606.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 15585.429 NA
1 28172.438 15585.429
2 34571.269 28172.438
3 33075.692 34571.269
4 42589.400 33075.692
5 21296.054 42589.400
6 24066.969 21296.054
7 27847.969 24066.969
8 30237.677 27847.969
9 27286.361 30237.677
10 40268.800 27286.361
11 25161.361 40268.800
12 31457.776 25161.361
13 41951.631 31457.776
14 26369.654 41951.631
15 44329.423 26369.654
16 45707.362 44329.423
17 32833.669 45707.362
18 34572.277 32833.669
19 45296.662 34572.277
20 27049.792 45296.662
21 56727.285 27049.792
22 45613.762 56727.285
23 33192.400 45613.762
24 38195.083 33192.400
25 36906.400 38195.083
26 23546.500 36906.400
27 23048.769 23546.500
28 5227.092 23048.769
29 14311.285 5227.092
30 2757.739 14311.285
31 3916.623 2757.739
32 -10197.361 3916.623
33 -6236.254 -10197.361
34 -23351.238 -6236.254
35 -16498.561 -23351.238
36 -15459.263 -16498.561
37 -46054.331 -15459.263
38 -24947.538 -46054.331
39 -32029.538 -24947.538
40 -41877.369 -32029.538
41 -19734.830 -41877.369
42 -23203.530 -19734.830
43 -53532.454 -23203.530
44 -22977.976 -53532.454
45 -48326.100 -22977.976
46 -47090.738 -48326.100
47 -36326.754 -47090.738
48 -55454.878 -36326.754
49 -60976.138 -55454.878
50 -59539.885 -60976.138
51 -68424.346 -59539.885
52 -51646.485 -68424.346
53 -48706.177 -51646.485
54 -38193.454 -48706.177
55 -23528.800 -38193.454
56 -24112.131 -23528.800
57 -29451.293 -24112.131
58 -15440.585 -29451.293
59 -5528.446 -15440.585
60 -14324.148 -5528.446
61 NA -14324.148
> dum1 <- dum[2:length(myerror),]
> dum1
lag(myerror, k = 1) myerror
[1,] 28172.438 15585.429
[2,] 34571.269 28172.438
[3,] 33075.692 34571.269
[4,] 42589.400 33075.692
[5,] 21296.054 42589.400
[6,] 24066.969 21296.054
[7,] 27847.969 24066.969
[8,] 30237.677 27847.969
[9,] 27286.361 30237.677
[10,] 40268.800 27286.361
[11,] 25161.361 40268.800
[12,] 31457.776 25161.361
[13,] 41951.631 31457.776
[14,] 26369.654 41951.631
[15,] 44329.423 26369.654
[16,] 45707.362 44329.423
[17,] 32833.669 45707.362
[18,] 34572.277 32833.669
[19,] 45296.662 34572.277
[20,] 27049.792 45296.662
[21,] 56727.285 27049.792
[22,] 45613.762 56727.285
[23,] 33192.400 45613.762
[24,] 38195.083 33192.400
[25,] 36906.400 38195.083
[26,] 23546.500 36906.400
[27,] 23048.769 23546.500
[28,] 5227.092 23048.769
[29,] 14311.285 5227.092
[30,] 2757.739 14311.285
[31,] 3916.623 2757.739
[32,] -10197.361 3916.623
[33,] -6236.254 -10197.361
[34,] -23351.238 -6236.254
[35,] -16498.561 -23351.238
[36,] -15459.263 -16498.561
[37,] -46054.331 -15459.263
[38,] -24947.538 -46054.331
[39,] -32029.538 -24947.538
[40,] -41877.369 -32029.538
[41,] -19734.830 -41877.369
[42,] -23203.530 -19734.830
[43,] -53532.454 -23203.530
[44,] -22977.976 -53532.454
[45,] -48326.100 -22977.976
[46,] -47090.738 -48326.100
[47,] -36326.754 -47090.738
[48,] -55454.878 -36326.754
[49,] -60976.138 -55454.878
[50,] -59539.885 -60976.138
[51,] -68424.346 -59539.885
[52,] -51646.485 -68424.346
[53,] -48706.177 -51646.485
[54,] -38193.454 -48706.177
[55,] -23528.800 -38193.454
[56,] -24112.131 -23528.800
[57,] -29451.293 -24112.131
[58,] -15440.585 -29451.293
[59,] -5528.446 -15440.585
[60,] -14324.148 -5528.446
> z <- as.data.frame(dum1)
> z
lag(myerror, k = 1) myerror
1 28172.438 15585.429
2 34571.269 28172.438
3 33075.692 34571.269
4 42589.400 33075.692
5 21296.054 42589.400
6 24066.969 21296.054
7 27847.969 24066.969
8 30237.677 27847.969
9 27286.361 30237.677
10 40268.800 27286.361
11 25161.361 40268.800
12 31457.776 25161.361
13 41951.631 31457.776
14 26369.654 41951.631
15 44329.423 26369.654
16 45707.362 44329.423
17 32833.669 45707.362
18 34572.277 32833.669
19 45296.662 34572.277
20 27049.792 45296.662
21 56727.285 27049.792
22 45613.762 56727.285
23 33192.400 45613.762
24 38195.083 33192.400
25 36906.400 38195.083
26 23546.500 36906.400
27 23048.769 23546.500
28 5227.092 23048.769
29 14311.285 5227.092
30 2757.739 14311.285
31 3916.623 2757.739
32 -10197.361 3916.623
33 -6236.254 -10197.361
34 -23351.238 -6236.254
35 -16498.561 -23351.238
36 -15459.263 -16498.561
37 -46054.331 -15459.263
38 -24947.538 -46054.331
39 -32029.538 -24947.538
40 -41877.369 -32029.538
41 -19734.830 -41877.369
42 -23203.530 -19734.830
43 -53532.454 -23203.530
44 -22977.976 -53532.454
45 -48326.100 -22977.976
46 -47090.738 -48326.100
47 -36326.754 -47090.738
48 -55454.878 -36326.754
49 -60976.138 -55454.878
50 -59539.885 -60976.138
51 -68424.346 -59539.885
52 -51646.485 -68424.346
53 -48706.177 -51646.485
54 -38193.454 -48706.177
55 -23528.800 -38193.454
56 -24112.131 -23528.800
57 -29451.293 -24112.131
58 -15440.585 -29451.293
59 -5528.446 -15440.585
60 -14324.148 -5528.446
> 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/7s3um1258623606.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/87g4u1258623606.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/9n6l01258623606.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/101z8w1258623606.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/1157rr1258623606.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/12159j1258623607.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/13y88i1258623607.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/1450qg1258623607.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/15c0601258623607.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/16mwrq1258623607.tab")
+ }
>
> system("convert tmp/1ymex1258623606.ps tmp/1ymex1258623606.png")
> system("convert tmp/2rfxb1258623606.ps tmp/2rfxb1258623606.png")
> system("convert tmp/33nob1258623606.ps tmp/33nob1258623606.png")
> system("convert tmp/45wgy1258623606.ps tmp/45wgy1258623606.png")
> system("convert tmp/52zbc1258623606.ps tmp/52zbc1258623606.png")
> system("convert tmp/63dl51258623606.ps tmp/63dl51258623606.png")
> system("convert tmp/7s3um1258623606.ps tmp/7s3um1258623606.png")
> system("convert tmp/87g4u1258623606.ps tmp/87g4u1258623606.png")
> system("convert tmp/9n6l01258623606.ps tmp/9n6l01258623606.png")
> system("convert tmp/101z8w1258623606.ps tmp/101z8w1258623606.png")
>
>
> proc.time()
user system elapsed
2.391 1.548 2.974