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 = '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 t
1 0 0 1
2 0 0 2
3 0 0 3
4 0 0 4
5 0 0 5
6 0 0 6
7 0 0 7
8 0 0 8
9 0 0 9
10 1 0 10
11 0 1 11
12 0 0 12
13 0 0 13
14 0 0 14
15 0 0 15
16 0 0 16
17 0 0 17
18 0 0 18
19 0 0 19
20 0 0 20
21 0 0 21
22 1 0 22
23 0 1 23
24 0 0 24
25 0 0 25
26 0 0 26
27 0 0 27
28 0 0 28
29 0 0 29
30 0 0 30
31 0 0 31
32 0 0 32
33 0 0 33
34 1 0 34
35 0 1 35
36 0 0 36
37 0 0 37
38 0 0 38
39 0 0 39
40 0 0 40
41 0 0 41
42 0 0 42
43 0 0 43
44 0 0 44
45 0 0 45
46 1 0 46
47 0 1 47
48 0 0 48
49 0 0 49
50 0 0 50
51 0 0 51
52 0 0 52
53 0 0 53
54 0 0 54
55 0 0 55
56 0 0 56
57 0 0 57
58 1 0 58
59 0 1 59
60 0 0 60
61 0 0 61
> k <- length(x[1,])
> df <- as.data.frame(x)
> (mylm <- lm(df))
Call:
lm(formula = df)
Coefficients:
(Intercept) Tot_ind_productie M1 M2
838838 -2239 25645 45056
M3 M4 M5 M6
31267 2946 -9777 -5649
M7 M8 M9 M10
-5301 7361 -9755 -20475
M11 t
3900 -1690
> (mysum <- summary(mylm))
Call:
lm(formula = df)
Residuals:
Min 1Q Median 3Q Max
-35288.4 -14386.1 -144.9 14481.4 37666.3
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 838838.2 44098.4 19.022 < 2e-16 ***
Tot_ind_productie -2239.1 484.4 -4.622 2.98e-05 ***
M1 25644.9 14160.6 1.811 0.0765 .
M2 45056.4 19007.2 2.370 0.0219 *
M3 31267.2 18958.9 1.649 0.1058
M4 2946.3 16970.2 0.174 0.8629
M5 -9776.7 15149.0 -0.645 0.5218
M6 -5649.1 15183.2 -0.372 0.7115
M7 -5301.2 15446.9 -0.343 0.7330
M8 7360.6 17738.5 0.415 0.6801
M9 -9755.2 16039.3 -0.608 0.5460
M10 -20474.8 15903.3 -1.287 0.2042
M11 3900.1 18513.8 0.211 0.8341
t -1690.1 161.7 -10.452 7.53e-14 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 21890 on 47 degrees of freedom
Multiple R-squared: 0.7923, Adjusted R-squared: 0.7348
F-statistic: 13.79 on 13 and 47 DF, p-value: 6.577e-12
> 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,] 4.494565e-02 8.989129e-02 9.550544e-01
[2,] 1.209485e-02 2.418969e-02 9.879052e-01
[3,] 3.016502e-03 6.033005e-03 9.969835e-01
[4,] 6.918966e-04 1.383793e-03 9.993081e-01
[5,] 2.512920e-04 5.025839e-04 9.997487e-01
[6,] 1.087094e-04 2.174189e-04 9.998913e-01
[7,] 2.531509e-05 5.063017e-05 9.999747e-01
[8,] 6.542737e-06 1.308547e-05 9.999935e-01
[9,] 1.284832e-06 2.569664e-06 9.999987e-01
[10,] 1.318652e-06 2.637304e-06 9.999987e-01
[11,] 9.578989e-05 1.915798e-04 9.999042e-01
[12,] 1.884868e-03 3.769735e-03 9.981151e-01
[13,] 2.634946e-02 5.269892e-02 9.736505e-01
[14,] 3.774218e-02 7.548436e-02 9.622578e-01
[15,] 5.726578e-02 1.145316e-01 9.427342e-01
[16,] 1.284191e-01 2.568382e-01 8.715809e-01
[17,] 1.522196e-01 3.044391e-01 8.477804e-01
[18,] 2.551414e-01 5.102827e-01 7.448586e-01
[19,] 3.730691e-01 7.461382e-01 6.269309e-01
[20,] 3.503538e-01 7.007076e-01 6.496462e-01
[21,] 3.114052e-01 6.228104e-01 6.885948e-01
[22,] 3.615608e-01 7.231216e-01 6.384392e-01
[23,] 5.162561e-01 9.674877e-01 4.837439e-01
[24,] 6.928936e-01 6.142128e-01 3.071064e-01
[25,] 9.013985e-01 1.972030e-01 9.860150e-02
[26,] 9.560807e-01 8.783856e-02 4.391928e-02
[27,] 9.857215e-01 2.855698e-02 1.427849e-02
[28,] 9.999824e-01 3.522565e-05 1.761282e-05
> postscript(file="/var/www/html/rcomp/tmp/143e41258625187.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/23jiv1258625187.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/35q241258625187.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/4jmr21258625187.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/5rlvz1258625187.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
-35288.3844 -12423.9991 -6742.1638 -8218.2572 2331.9766 -20368.0835
7 8 9 10 11 12
-17601.0649 -13820.0649 -10596.4603 -14354.3956 -89.9380 -16682.0247
13 14 15 16 17 18
650.3261 21460.5246 4479.6271 24270.8522 25867.0067 12307.4926
19 20 21 22 23 24
14120.1379 26013.5369 6087.1837 37666.2727 25672.0921 12506.4581
25 26 27 28 29 30
28603.5288 37007.7146 24002.4156 24069.7085 5608.9703 15429.6417
31 32 33 34 35 36
3482.5274 5557.1399 -9827.1737 -5133.4836 -22505.6514 -15520.4862
37 38 39 40 41 42
-4049.2415 -26120.4548 -3061.4085 -10649.9814 -21604.4791 2579.9378
43 44 45 46 47 48
-736.7903 -33364.7750 -144.9453 -26981.6133 -26217.7543 -14386.0699
49 50 51 52 53 54
-24621.4505 -19923.7853 -18678.4704 -29472.3221 -12203.4745 -9948.9886
55 56 57 58 59 60
735.1899 15614.1631 14481.3955 8803.2198 23141.2515 34082.1227
61
34705.2215
> postscript(file="/var/www/html/rcomp/tmp/6t5w71258625187.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 -35288.3844 NA
1 -12423.9991 -35288.3844
2 -6742.1638 -12423.9991
3 -8218.2572 -6742.1638
4 2331.9766 -8218.2572
5 -20368.0835 2331.9766
6 -17601.0649 -20368.0835
7 -13820.0649 -17601.0649
8 -10596.4603 -13820.0649
9 -14354.3956 -10596.4603
10 -89.9380 -14354.3956
11 -16682.0247 -89.9380
12 650.3261 -16682.0247
13 21460.5246 650.3261
14 4479.6271 21460.5246
15 24270.8522 4479.6271
16 25867.0067 24270.8522
17 12307.4926 25867.0067
18 14120.1379 12307.4926
19 26013.5369 14120.1379
20 6087.1837 26013.5369
21 37666.2727 6087.1837
22 25672.0921 37666.2727
23 12506.4581 25672.0921
24 28603.5288 12506.4581
25 37007.7146 28603.5288
26 24002.4156 37007.7146
27 24069.7085 24002.4156
28 5608.9703 24069.7085
29 15429.6417 5608.9703
30 3482.5274 15429.6417
31 5557.1399 3482.5274
32 -9827.1737 5557.1399
33 -5133.4836 -9827.1737
34 -22505.6514 -5133.4836
35 -15520.4862 -22505.6514
36 -4049.2415 -15520.4862
37 -26120.4548 -4049.2415
38 -3061.4085 -26120.4548
39 -10649.9814 -3061.4085
40 -21604.4791 -10649.9814
41 2579.9378 -21604.4791
42 -736.7903 2579.9378
43 -33364.7750 -736.7903
44 -144.9453 -33364.7750
45 -26981.6133 -144.9453
46 -26217.7543 -26981.6133
47 -14386.0699 -26217.7543
48 -24621.4505 -14386.0699
49 -19923.7853 -24621.4505
50 -18678.4704 -19923.7853
51 -29472.3221 -18678.4704
52 -12203.4745 -29472.3221
53 -9948.9886 -12203.4745
54 735.1899 -9948.9886
55 15614.1631 735.1899
56 14481.3955 15614.1631
57 8803.2198 14481.3955
58 23141.2515 8803.2198
59 34082.1227 23141.2515
60 34705.2215 34082.1227
61 NA 34705.2215
> dum1 <- dum[2:length(myerror),]
> dum1
lag(myerror, k = 1) myerror
[1,] -12423.9991 -35288.3844
[2,] -6742.1638 -12423.9991
[3,] -8218.2572 -6742.1638
[4,] 2331.9766 -8218.2572
[5,] -20368.0835 2331.9766
[6,] -17601.0649 -20368.0835
[7,] -13820.0649 -17601.0649
[8,] -10596.4603 -13820.0649
[9,] -14354.3956 -10596.4603
[10,] -89.9380 -14354.3956
[11,] -16682.0247 -89.9380
[12,] 650.3261 -16682.0247
[13,] 21460.5246 650.3261
[14,] 4479.6271 21460.5246
[15,] 24270.8522 4479.6271
[16,] 25867.0067 24270.8522
[17,] 12307.4926 25867.0067
[18,] 14120.1379 12307.4926
[19,] 26013.5369 14120.1379
[20,] 6087.1837 26013.5369
[21,] 37666.2727 6087.1837
[22,] 25672.0921 37666.2727
[23,] 12506.4581 25672.0921
[24,] 28603.5288 12506.4581
[25,] 37007.7146 28603.5288
[26,] 24002.4156 37007.7146
[27,] 24069.7085 24002.4156
[28,] 5608.9703 24069.7085
[29,] 15429.6417 5608.9703
[30,] 3482.5274 15429.6417
[31,] 5557.1399 3482.5274
[32,] -9827.1737 5557.1399
[33,] -5133.4836 -9827.1737
[34,] -22505.6514 -5133.4836
[35,] -15520.4862 -22505.6514
[36,] -4049.2415 -15520.4862
[37,] -26120.4548 -4049.2415
[38,] -3061.4085 -26120.4548
[39,] -10649.9814 -3061.4085
[40,] -21604.4791 -10649.9814
[41,] 2579.9378 -21604.4791
[42,] -736.7903 2579.9378
[43,] -33364.7750 -736.7903
[44,] -144.9453 -33364.7750
[45,] -26981.6133 -144.9453
[46,] -26217.7543 -26981.6133
[47,] -14386.0699 -26217.7543
[48,] -24621.4505 -14386.0699
[49,] -19923.7853 -24621.4505
[50,] -18678.4704 -19923.7853
[51,] -29472.3221 -18678.4704
[52,] -12203.4745 -29472.3221
[53,] -9948.9886 -12203.4745
[54,] 735.1899 -9948.9886
[55,] 15614.1631 735.1899
[56,] 14481.3955 15614.1631
[57,] 8803.2198 14481.3955
[58,] 23141.2515 8803.2198
[59,] 34082.1227 23141.2515
[60,] 34705.2215 34082.1227
> z <- as.data.frame(dum1)
> z
lag(myerror, k = 1) myerror
1 -12423.9991 -35288.3844
2 -6742.1638 -12423.9991
3 -8218.2572 -6742.1638
4 2331.9766 -8218.2572
5 -20368.0835 2331.9766
6 -17601.0649 -20368.0835
7 -13820.0649 -17601.0649
8 -10596.4603 -13820.0649
9 -14354.3956 -10596.4603
10 -89.9380 -14354.3956
11 -16682.0247 -89.9380
12 650.3261 -16682.0247
13 21460.5246 650.3261
14 4479.6271 21460.5246
15 24270.8522 4479.6271
16 25867.0067 24270.8522
17 12307.4926 25867.0067
18 14120.1379 12307.4926
19 26013.5369 14120.1379
20 6087.1837 26013.5369
21 37666.2727 6087.1837
22 25672.0921 37666.2727
23 12506.4581 25672.0921
24 28603.5288 12506.4581
25 37007.7146 28603.5288
26 24002.4156 37007.7146
27 24069.7085 24002.4156
28 5608.9703 24069.7085
29 15429.6417 5608.9703
30 3482.5274 15429.6417
31 5557.1399 3482.5274
32 -9827.1737 5557.1399
33 -5133.4836 -9827.1737
34 -22505.6514 -5133.4836
35 -15520.4862 -22505.6514
36 -4049.2415 -15520.4862
37 -26120.4548 -4049.2415
38 -3061.4085 -26120.4548
39 -10649.9814 -3061.4085
40 -21604.4791 -10649.9814
41 2579.9378 -21604.4791
42 -736.7903 2579.9378
43 -33364.7750 -736.7903
44 -144.9453 -33364.7750
45 -26981.6133 -144.9453
46 -26217.7543 -26981.6133
47 -14386.0699 -26217.7543
48 -24621.4505 -14386.0699
49 -19923.7853 -24621.4505
50 -18678.4704 -19923.7853
51 -29472.3221 -18678.4704
52 -12203.4745 -29472.3221
53 -9948.9886 -12203.4745
54 735.1899 -9948.9886
55 15614.1631 735.1899
56 14481.3955 15614.1631
57 8803.2198 14481.3955
58 23141.2515 8803.2198
59 34082.1227 23141.2515
60 34705.2215 34082.1227
> 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/7a0tq1258625187.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/88ez21258625187.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/97psm1258625187.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/10encb1258625187.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/11fx5w1258625187.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/12smcy1258625187.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/13w0vx1258625187.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/14v8px1258625187.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/150dqw1258625188.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/16sp2m1258625188.tab")
+ }
>
> system("convert tmp/143e41258625187.ps tmp/143e41258625187.png")
> system("convert tmp/23jiv1258625187.ps tmp/23jiv1258625187.png")
> system("convert tmp/35q241258625187.ps tmp/35q241258625187.png")
> system("convert tmp/4jmr21258625187.ps tmp/4jmr21258625187.png")
> system("convert tmp/5rlvz1258625187.ps tmp/5rlvz1258625187.png")
> system("convert tmp/6t5w71258625187.ps tmp/6t5w71258625187.png")
> system("convert tmp/7a0tq1258625187.ps tmp/7a0tq1258625187.png")
> system("convert tmp/88ez21258625187.ps tmp/88ez21258625187.png")
> system("convert tmp/97psm1258625187.ps tmp/97psm1258625187.png")
> system("convert tmp/10encb1258625187.ps tmp/10encb1258625187.png")
>
>
> proc.time()
user system elapsed
2.325 1.491 3.656