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(7291,4071,6820,4351,8031,4871,7862,4649,7357,4922,7213,4879,7079,4853,7012,4545,7319,4733,8148,5191,7599,4983,6908,4593,7878,4656,7407,4513,7911,4857,7323,4681,7179,4897,6758,4547,6934,4692,6696,4390,7688,5341,8296,5415,7697,4890,7907,5120,7592,4422,7710,4797,9011,5689,8225,5171,7733,4265,8062,5215,7859,4874,8221,4590,8330,4994,8868,4988,9053,5110,8811,5141,8120,4395,7953,4523,8878,5306,8601,5365,8361,5496,9116,5647,9310,5443,9891,5546,10147,5912,10317,5665,10682,5963,10276,5861,10614,5366,9413,5619,11068,6721,9772,6054,10350,6619,10541,6856,10049,6193,10714,6317,10759,6618,11684,6585,11462,6852,10485,6586,11056,6154,10184,6193,11082,7606,10554,6588,11315,7143,10847,7629,11104,7041,11026,7146,11073,7200,12073,7739,12328,7953,11172,7082),dim=c(2,72),dimnames=list(c('UitvEU','Uitvniet-EU'),1:72))
> y <- array(NA,dim=c(2,72),dimnames=list(c('UitvEU','Uitvniet-EU'),1:72))
> 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
UitvEU Uitvniet-EU M1 M2 M3 M4 M5 M6 M7 M8 M9 M10 M11
1 7291 4071 1 0 0 0 0 0 0 0 0 0 0
2 6820 4351 0 1 0 0 0 0 0 0 0 0 0
3 8031 4871 0 0 1 0 0 0 0 0 0 0 0
4 7862 4649 0 0 0 1 0 0 0 0 0 0 0
5 7357 4922 0 0 0 0 1 0 0 0 0 0 0
6 7213 4879 0 0 0 0 0 1 0 0 0 0 0
7 7079 4853 0 0 0 0 0 0 1 0 0 0 0
8 7012 4545 0 0 0 0 0 0 0 1 0 0 0
9 7319 4733 0 0 0 0 0 0 0 0 1 0 0
10 8148 5191 0 0 0 0 0 0 0 0 0 1 0
11 7599 4983 0 0 0 0 0 0 0 0 0 0 1
12 6908 4593 0 0 0 0 0 0 0 0 0 0 0
13 7878 4656 1 0 0 0 0 0 0 0 0 0 0
14 7407 4513 0 1 0 0 0 0 0 0 0 0 0
15 7911 4857 0 0 1 0 0 0 0 0 0 0 0
16 7323 4681 0 0 0 1 0 0 0 0 0 0 0
17 7179 4897 0 0 0 0 1 0 0 0 0 0 0
18 6758 4547 0 0 0 0 0 1 0 0 0 0 0
19 6934 4692 0 0 0 0 0 0 1 0 0 0 0
20 6696 4390 0 0 0 0 0 0 0 1 0 0 0
21 7688 5341 0 0 0 0 0 0 0 0 1 0 0
22 8296 5415 0 0 0 0 0 0 0 0 0 1 0
23 7697 4890 0 0 0 0 0 0 0 0 0 0 1
24 7907 5120 0 0 0 0 0 0 0 0 0 0 0
25 7592 4422 1 0 0 0 0 0 0 0 0 0 0
26 7710 4797 0 1 0 0 0 0 0 0 0 0 0
27 9011 5689 0 0 1 0 0 0 0 0 0 0 0
28 8225 5171 0 0 0 1 0 0 0 0 0 0 0
29 7733 4265 0 0 0 0 1 0 0 0 0 0 0
30 8062 5215 0 0 0 0 0 1 0 0 0 0 0
31 7859 4874 0 0 0 0 0 0 1 0 0 0 0
32 8221 4590 0 0 0 0 0 0 0 1 0 0 0
33 8330 4994 0 0 0 0 0 0 0 0 1 0 0
34 8868 4988 0 0 0 0 0 0 0 0 0 1 0
35 9053 5110 0 0 0 0 0 0 0 0 0 0 1
36 8811 5141 0 0 0 0 0 0 0 0 0 0 0
37 8120 4395 1 0 0 0 0 0 0 0 0 0 0
38 7953 4523 0 1 0 0 0 0 0 0 0 0 0
39 8878 5306 0 0 1 0 0 0 0 0 0 0 0
40 8601 5365 0 0 0 1 0 0 0 0 0 0 0
41 8361 5496 0 0 0 0 1 0 0 0 0 0 0
42 9116 5647 0 0 0 0 0 1 0 0 0 0 0
43 9310 5443 0 0 0 0 0 0 1 0 0 0 0
44 9891 5546 0 0 0 0 0 0 0 1 0 0 0
45 10147 5912 0 0 0 0 0 0 0 0 1 0 0
46 10317 5665 0 0 0 0 0 0 0 0 0 1 0
47 10682 5963 0 0 0 0 0 0 0 0 0 0 1
48 10276 5861 0 0 0 0 0 0 0 0 0 0 0
49 10614 5366 1 0 0 0 0 0 0 0 0 0 0
50 9413 5619 0 1 0 0 0 0 0 0 0 0 0
51 11068 6721 0 0 1 0 0 0 0 0 0 0 0
52 9772 6054 0 0 0 1 0 0 0 0 0 0 0
53 10350 6619 0 0 0 0 1 0 0 0 0 0 0
54 10541 6856 0 0 0 0 0 1 0 0 0 0 0
55 10049 6193 0 0 0 0 0 0 1 0 0 0 0
56 10714 6317 0 0 0 0 0 0 0 1 0 0 0
57 10759 6618 0 0 0 0 0 0 0 0 1 0 0
58 11684 6585 0 0 0 0 0 0 0 0 0 1 0
59 11462 6852 0 0 0 0 0 0 0 0 0 0 1
60 10485 6586 0 0 0 0 0 0 0 0 0 0 0
61 11056 6154 1 0 0 0 0 0 0 0 0 0 0
62 10184 6193 0 1 0 0 0 0 0 0 0 0 0
63 11082 7606 0 0 1 0 0 0 0 0 0 0 0
64 10554 6588 0 0 0 1 0 0 0 0 0 0 0
65 11315 7143 0 0 0 0 1 0 0 0 0 0 0
66 10847 7629 0 0 0 0 0 1 0 0 0 0 0
67 11104 7041 0 0 0 0 0 0 1 0 0 0 0
68 11026 7146 0 0 0 0 0 0 0 1 0 0 0
69 11073 7200 0 0 0 0 0 0 0 0 1 0 0
70 12073 7739 0 0 0 0 0 0 0 0 0 1 0
71 12328 7953 0 0 0 0 0 0 0 0 0 0 1
72 11172 7082 0 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) `Uitvniet-EU` M1 M2 M3
451.553 1.537 861.295 111.867 -100.540
M4 M5 M6 M7 M8
-56.660 -277.315 -603.577 -207.628 140.512
M9 M10 M11
-146.815 330.416 193.211
> (mysum <- summary(mylm))
Call:
lm(formula = df)
Residuals:
Min 1Q Median 3Q Max
-960.10 -439.50 17.28 354.59 1053.14
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 451.55324 460.57506 0.980 0.33089
`Uitvniet-EU` 1.53709 0.07021 21.892 < 2e-16 ***
M1 861.29476 323.02543 2.666 0.00988 **
M2 111.86717 321.10221 0.348 0.72879
M3 -100.53956 317.06776 -0.317 0.75229
M4 -56.66015 317.73021 -0.178 0.85908
M5 -277.31531 317.20568 -0.874 0.38553
M6 -603.57736 317.00452 -1.904 0.06179 .
M7 -207.62806 317.32928 -0.654 0.51546
M8 140.51247 317.70934 0.442 0.65991
M9 -146.81522 317.00887 -0.463 0.64498
M10 330.41583 317.28259 1.041 0.30194
M11 193.21071 317.37568 0.609 0.54501
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 549 on 59 degrees of freedom
Multiple R-squared: 0.9, Adjusted R-squared: 0.8796
F-statistic: 44.23 on 12 and 59 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.1525420037 0.305084007 0.847457996
[2,] 0.0712121258 0.142424252 0.928787874
[3,] 0.0278196024 0.055639205 0.972180398
[4,] 0.0114777863 0.022955573 0.988522214
[5,] 0.0056563236 0.011312647 0.994343676
[6,] 0.0046380479 0.009276096 0.995361952
[7,] 0.0030731960 0.006146392 0.996926804
[8,] 0.0023859593 0.004771919 0.997614041
[9,] 0.0044155130 0.008831026 0.995584487
[10,] 0.0046004838 0.009200968 0.995399516
[11,] 0.0033668741 0.006733748 0.996633126
[12,] 0.0015704090 0.003140818 0.998429591
[13,] 0.0007559605 0.001511921 0.999244039
[14,] 0.0841499452 0.168299890 0.915850055
[15,] 0.1077699940 0.215539988 0.892230006
[16,] 0.2004068530 0.400813706 0.799593147
[17,] 0.4706789987 0.941357997 0.529321001
[18,] 0.5688263839 0.862347232 0.431173616
[19,] 0.6886642471 0.622671506 0.311335753
[20,] 0.8094397615 0.381120477 0.190560239
[21,] 0.8538634800 0.292273040 0.146136520
[22,] 0.9552878247 0.089424351 0.044712175
[23,] 0.9574027848 0.085194430 0.042597215
[24,] 0.9474790035 0.105041993 0.052520996
[25,] 0.9533576054 0.093284789 0.046642395
[26,] 0.9973025898 0.005394820 0.002697410
[27,] 0.9970913462 0.005817308 0.002908654
[28,] 0.9979466875 0.004106625 0.002053312
[29,] 0.9970764859 0.005847028 0.002923514
[30,] 0.9952735676 0.009452865 0.004726432
[31,] 0.9965854317 0.006829137 0.003414568
[32,] 0.9946158222 0.010768356 0.005384178
[33,] 0.9897282892 0.020543422 0.010271711
[34,] 0.9826579157 0.034684169 0.017342084
[35,] 0.9747754917 0.050449017 0.025224508
[36,] 0.9793177366 0.041364527 0.020682263
[37,] 0.9707534690 0.058493062 0.029246531
[38,] 0.9809995907 0.038000819 0.019000409
[39,] 0.9604611402 0.079077720 0.039538860
[40,] 0.9671700697 0.065659861 0.032829930
[41,] 0.9164246753 0.167150649 0.083575325
> postscript(file="/var/www/html/rcomp/tmp/1h5rm1258562722.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/2ap6n1258562722.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/3ujc11258562722.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/431o31258562722.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/520bn1258562722.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 = 72
Frequency = 1
1 2 3 4 5 6
-279.331339 -431.288251 192.832963 321.186978 -382.782751 -134.425945
7 8 9 10 11 12
-624.410967 -566.128540 -260.773301 -612.990438 -705.071114 -603.396267
13 14 15 16 17 18
-591.527544 -93.296431 94.352189 -266.999823 -522.355563 -79.112885
19 20 21 22 23 24
-521.939875 -643.879973 -826.322519 -809.298045 -464.121974 -414.441396
25 26 27 28 29 30
-517.849062 -226.829290 -84.504637 -118.172713 1003.083756 198.112645
31 32 33 34 35 36
123.310195 573.702521 349.046854 419.038331 553.718770 457.279766
37 38 39 40 41 42
51.652301 437.332694 371.199887 -40.367694 -261.070994 588.090832
43 44 45 46 47 48
699.707390 774.246842 755.000501 827.430073 871.583106 815.576744
49 50 51 52 53 54
1053.140309 212.684760 386.221031 71.578997 1.779709 154.752007
55 56 57 58 59 60
285.891742 412.152356 281.816704 780.309545 285.112291 -89.811716
61 62 63 64 65 66
283.915335 101.396518 -960.101433 32.774256 161.345843 -727.416654
67 68 69 70 71 72
37.441516 -550.093207 -298.768239 -604.489466 -541.221080 -165.207131
> postscript(file="/var/www/html/rcomp/tmp/62cob1258562722.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 = 72
Frequency = 1
lag(myerror, k = 1) myerror
0 -279.331339 NA
1 -431.288251 -279.331339
2 192.832963 -431.288251
3 321.186978 192.832963
4 -382.782751 321.186978
5 -134.425945 -382.782751
6 -624.410967 -134.425945
7 -566.128540 -624.410967
8 -260.773301 -566.128540
9 -612.990438 -260.773301
10 -705.071114 -612.990438
11 -603.396267 -705.071114
12 -591.527544 -603.396267
13 -93.296431 -591.527544
14 94.352189 -93.296431
15 -266.999823 94.352189
16 -522.355563 -266.999823
17 -79.112885 -522.355563
18 -521.939875 -79.112885
19 -643.879973 -521.939875
20 -826.322519 -643.879973
21 -809.298045 -826.322519
22 -464.121974 -809.298045
23 -414.441396 -464.121974
24 -517.849062 -414.441396
25 -226.829290 -517.849062
26 -84.504637 -226.829290
27 -118.172713 -84.504637
28 1003.083756 -118.172713
29 198.112645 1003.083756
30 123.310195 198.112645
31 573.702521 123.310195
32 349.046854 573.702521
33 419.038331 349.046854
34 553.718770 419.038331
35 457.279766 553.718770
36 51.652301 457.279766
37 437.332694 51.652301
38 371.199887 437.332694
39 -40.367694 371.199887
40 -261.070994 -40.367694
41 588.090832 -261.070994
42 699.707390 588.090832
43 774.246842 699.707390
44 755.000501 774.246842
45 827.430073 755.000501
46 871.583106 827.430073
47 815.576744 871.583106
48 1053.140309 815.576744
49 212.684760 1053.140309
50 386.221031 212.684760
51 71.578997 386.221031
52 1.779709 71.578997
53 154.752007 1.779709
54 285.891742 154.752007
55 412.152356 285.891742
56 281.816704 412.152356
57 780.309545 281.816704
58 285.112291 780.309545
59 -89.811716 285.112291
60 283.915335 -89.811716
61 101.396518 283.915335
62 -960.101433 101.396518
63 32.774256 -960.101433
64 161.345843 32.774256
65 -727.416654 161.345843
66 37.441516 -727.416654
67 -550.093207 37.441516
68 -298.768239 -550.093207
69 -604.489466 -298.768239
70 -541.221080 -604.489466
71 -165.207131 -541.221080
72 NA -165.207131
> dum1 <- dum[2:length(myerror),]
> dum1
lag(myerror, k = 1) myerror
[1,] -431.288251 -279.331339
[2,] 192.832963 -431.288251
[3,] 321.186978 192.832963
[4,] -382.782751 321.186978
[5,] -134.425945 -382.782751
[6,] -624.410967 -134.425945
[7,] -566.128540 -624.410967
[8,] -260.773301 -566.128540
[9,] -612.990438 -260.773301
[10,] -705.071114 -612.990438
[11,] -603.396267 -705.071114
[12,] -591.527544 -603.396267
[13,] -93.296431 -591.527544
[14,] 94.352189 -93.296431
[15,] -266.999823 94.352189
[16,] -522.355563 -266.999823
[17,] -79.112885 -522.355563
[18,] -521.939875 -79.112885
[19,] -643.879973 -521.939875
[20,] -826.322519 -643.879973
[21,] -809.298045 -826.322519
[22,] -464.121974 -809.298045
[23,] -414.441396 -464.121974
[24,] -517.849062 -414.441396
[25,] -226.829290 -517.849062
[26,] -84.504637 -226.829290
[27,] -118.172713 -84.504637
[28,] 1003.083756 -118.172713
[29,] 198.112645 1003.083756
[30,] 123.310195 198.112645
[31,] 573.702521 123.310195
[32,] 349.046854 573.702521
[33,] 419.038331 349.046854
[34,] 553.718770 419.038331
[35,] 457.279766 553.718770
[36,] 51.652301 457.279766
[37,] 437.332694 51.652301
[38,] 371.199887 437.332694
[39,] -40.367694 371.199887
[40,] -261.070994 -40.367694
[41,] 588.090832 -261.070994
[42,] 699.707390 588.090832
[43,] 774.246842 699.707390
[44,] 755.000501 774.246842
[45,] 827.430073 755.000501
[46,] 871.583106 827.430073
[47,] 815.576744 871.583106
[48,] 1053.140309 815.576744
[49,] 212.684760 1053.140309
[50,] 386.221031 212.684760
[51,] 71.578997 386.221031
[52,] 1.779709 71.578997
[53,] 154.752007 1.779709
[54,] 285.891742 154.752007
[55,] 412.152356 285.891742
[56,] 281.816704 412.152356
[57,] 780.309545 281.816704
[58,] 285.112291 780.309545
[59,] -89.811716 285.112291
[60,] 283.915335 -89.811716
[61,] 101.396518 283.915335
[62,] -960.101433 101.396518
[63,] 32.774256 -960.101433
[64,] 161.345843 32.774256
[65,] -727.416654 161.345843
[66,] 37.441516 -727.416654
[67,] -550.093207 37.441516
[68,] -298.768239 -550.093207
[69,] -604.489466 -298.768239
[70,] -541.221080 -604.489466
[71,] -165.207131 -541.221080
> z <- as.data.frame(dum1)
> z
lag(myerror, k = 1) myerror
1 -431.288251 -279.331339
2 192.832963 -431.288251
3 321.186978 192.832963
4 -382.782751 321.186978
5 -134.425945 -382.782751
6 -624.410967 -134.425945
7 -566.128540 -624.410967
8 -260.773301 -566.128540
9 -612.990438 -260.773301
10 -705.071114 -612.990438
11 -603.396267 -705.071114
12 -591.527544 -603.396267
13 -93.296431 -591.527544
14 94.352189 -93.296431
15 -266.999823 94.352189
16 -522.355563 -266.999823
17 -79.112885 -522.355563
18 -521.939875 -79.112885
19 -643.879973 -521.939875
20 -826.322519 -643.879973
21 -809.298045 -826.322519
22 -464.121974 -809.298045
23 -414.441396 -464.121974
24 -517.849062 -414.441396
25 -226.829290 -517.849062
26 -84.504637 -226.829290
27 -118.172713 -84.504637
28 1003.083756 -118.172713
29 198.112645 1003.083756
30 123.310195 198.112645
31 573.702521 123.310195
32 349.046854 573.702521
33 419.038331 349.046854
34 553.718770 419.038331
35 457.279766 553.718770
36 51.652301 457.279766
37 437.332694 51.652301
38 371.199887 437.332694
39 -40.367694 371.199887
40 -261.070994 -40.367694
41 588.090832 -261.070994
42 699.707390 588.090832
43 774.246842 699.707390
44 755.000501 774.246842
45 827.430073 755.000501
46 871.583106 827.430073
47 815.576744 871.583106
48 1053.140309 815.576744
49 212.684760 1053.140309
50 386.221031 212.684760
51 71.578997 386.221031
52 1.779709 71.578997
53 154.752007 1.779709
54 285.891742 154.752007
55 412.152356 285.891742
56 281.816704 412.152356
57 780.309545 281.816704
58 285.112291 780.309545
59 -89.811716 285.112291
60 283.915335 -89.811716
61 101.396518 283.915335
62 -960.101433 101.396518
63 32.774256 -960.101433
64 161.345843 32.774256
65 -727.416654 161.345843
66 37.441516 -727.416654
67 -550.093207 37.441516
68 -298.768239 -550.093207
69 -604.489466 -298.768239
70 -541.221080 -604.489466
71 -165.207131 -541.221080
> 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/73laa1258562722.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/8nmq51258562722.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/9ecxz1258562722.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/106nei1258562722.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/11f4rm1258562722.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/12tazu1258562722.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/13itlh1258562722.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/14lwok1258562722.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/15z3051258562722.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/16ul921258562722.tab")
+ }
>
> system("convert tmp/1h5rm1258562722.ps tmp/1h5rm1258562722.png")
> system("convert tmp/2ap6n1258562722.ps tmp/2ap6n1258562722.png")
> system("convert tmp/3ujc11258562722.ps tmp/3ujc11258562722.png")
> system("convert tmp/431o31258562722.ps tmp/431o31258562722.png")
> system("convert tmp/520bn1258562722.ps tmp/520bn1258562722.png")
> system("convert tmp/62cob1258562722.ps tmp/62cob1258562722.png")
> system("convert tmp/73laa1258562722.ps tmp/73laa1258562722.png")
> system("convert tmp/8nmq51258562722.ps tmp/8nmq51258562722.png")
> system("convert tmp/9ecxz1258562722.ps tmp/9ecxz1258562722.png")
> system("convert tmp/106nei1258562722.ps tmp/106nei1258562722.png")
>
>
> proc.time()
user system elapsed
2.562 1.558 3.047