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 = '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 t
1 7291 4071 1 0 0 0 0 0 0 0 0 0 0 1
2 6820 4351 0 1 0 0 0 0 0 0 0 0 0 2
3 8031 4871 0 0 1 0 0 0 0 0 0 0 0 3
4 7862 4649 0 0 0 1 0 0 0 0 0 0 0 4
5 7357 4922 0 0 0 0 1 0 0 0 0 0 0 5
6 7213 4879 0 0 0 0 0 1 0 0 0 0 0 6
7 7079 4853 0 0 0 0 0 0 1 0 0 0 0 7
8 7012 4545 0 0 0 0 0 0 0 1 0 0 0 8
9 7319 4733 0 0 0 0 0 0 0 0 1 0 0 9
10 8148 5191 0 0 0 0 0 0 0 0 0 1 0 10
11 7599 4983 0 0 0 0 0 0 0 0 0 0 1 11
12 6908 4593 0 0 0 0 0 0 0 0 0 0 0 12
13 7878 4656 1 0 0 0 0 0 0 0 0 0 0 13
14 7407 4513 0 1 0 0 0 0 0 0 0 0 0 14
15 7911 4857 0 0 1 0 0 0 0 0 0 0 0 15
16 7323 4681 0 0 0 1 0 0 0 0 0 0 0 16
17 7179 4897 0 0 0 0 1 0 0 0 0 0 0 17
18 6758 4547 0 0 0 0 0 1 0 0 0 0 0 18
19 6934 4692 0 0 0 0 0 0 1 0 0 0 0 19
20 6696 4390 0 0 0 0 0 0 0 1 0 0 0 20
21 7688 5341 0 0 0 0 0 0 0 0 1 0 0 21
22 8296 5415 0 0 0 0 0 0 0 0 0 1 0 22
23 7697 4890 0 0 0 0 0 0 0 0 0 0 1 23
24 7907 5120 0 0 0 0 0 0 0 0 0 0 0 24
25 7592 4422 1 0 0 0 0 0 0 0 0 0 0 25
26 7710 4797 0 1 0 0 0 0 0 0 0 0 0 26
27 9011 5689 0 0 1 0 0 0 0 0 0 0 0 27
28 8225 5171 0 0 0 1 0 0 0 0 0 0 0 28
29 7733 4265 0 0 0 0 1 0 0 0 0 0 0 29
30 8062 5215 0 0 0 0 0 1 0 0 0 0 0 30
31 7859 4874 0 0 0 0 0 0 1 0 0 0 0 31
32 8221 4590 0 0 0 0 0 0 0 1 0 0 0 32
33 8330 4994 0 0 0 0 0 0 0 0 1 0 0 33
34 8868 4988 0 0 0 0 0 0 0 0 0 1 0 34
35 9053 5110 0 0 0 0 0 0 0 0 0 0 1 35
36 8811 5141 0 0 0 0 0 0 0 0 0 0 0 36
37 8120 4395 1 0 0 0 0 0 0 0 0 0 0 37
38 7953 4523 0 1 0 0 0 0 0 0 0 0 0 38
39 8878 5306 0 0 1 0 0 0 0 0 0 0 0 39
40 8601 5365 0 0 0 1 0 0 0 0 0 0 0 40
41 8361 5496 0 0 0 0 1 0 0 0 0 0 0 41
42 9116 5647 0 0 0 0 0 1 0 0 0 0 0 42
43 9310 5443 0 0 0 0 0 0 1 0 0 0 0 43
44 9891 5546 0 0 0 0 0 0 0 1 0 0 0 44
45 10147 5912 0 0 0 0 0 0 0 0 1 0 0 45
46 10317 5665 0 0 0 0 0 0 0 0 0 1 0 46
47 10682 5963 0 0 0 0 0 0 0 0 0 0 1 47
48 10276 5861 0 0 0 0 0 0 0 0 0 0 0 48
49 10614 5366 1 0 0 0 0 0 0 0 0 0 0 49
50 9413 5619 0 1 0 0 0 0 0 0 0 0 0 50
51 11068 6721 0 0 1 0 0 0 0 0 0 0 0 51
52 9772 6054 0 0 0 1 0 0 0 0 0 0 0 52
53 10350 6619 0 0 0 0 1 0 0 0 0 0 0 53
54 10541 6856 0 0 0 0 0 1 0 0 0 0 0 54
55 10049 6193 0 0 0 0 0 0 1 0 0 0 0 55
56 10714 6317 0 0 0 0 0 0 0 1 0 0 0 56
57 10759 6618 0 0 0 0 0 0 0 0 1 0 0 57
58 11684 6585 0 0 0 0 0 0 0 0 0 1 0 58
59 11462 6852 0 0 0 0 0 0 0 0 0 0 1 59
60 10485 6586 0 0 0 0 0 0 0 0 0 0 0 60
61 11056 6154 1 0 0 0 0 0 0 0 0 0 0 61
62 10184 6193 0 1 0 0 0 0 0 0 0 0 0 62
63 11082 7606 0 0 1 0 0 0 0 0 0 0 0 63
64 10554 6588 0 0 0 1 0 0 0 0 0 0 0 64
65 11315 7143 0 0 0 0 1 0 0 0 0 0 0 65
66 10847 7629 0 0 0 0 0 1 0 0 0 0 0 66
67 11104 7041 0 0 0 0 0 0 1 0 0 0 0 67
68 11026 7146 0 0 0 0 0 0 0 1 0 0 0 68
69 11073 7200 0 0 0 0 0 0 0 0 1 0 0 69
70 12073 7739 0 0 0 0 0 0 0 0 0 1 0 70
71 12328 7953 0 0 0 0 0 0 0 0 0 0 1 71
72 11172 7082 0 0 0 0 0 0 0 0 0 0 0 72
> k <- length(x[1,])
> df <- as.data.frame(x)
> (mylm <- lm(df))
Call:
lm(formula = df)
Coefficients:
(Intercept) `Uitvniet-EU` M1 M2 M3
3229.181 0.775 602.002 -66.893 324.782
M4 M5 M6 M7 M8
7.948 -144.623 -326.973 -181.870 57.044
M9 M10 M11 t
19.430 558.522 404.810 37.845
> (mysum <- summary(mylm))
Call:
lm(formula = df)
Residuals:
Min 1Q Median 3Q Max
-750.88 -279.91 -27.54 246.12 870.49
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 3229.181 539.659 5.984 1.45e-07 ***
`Uitvniet-EU` 0.775 0.125 6.201 6.34e-08 ***
M1 602.002 247.086 2.436 0.0179 *
M2 -66.893 244.062 -0.274 0.7850
M3 324.782 247.752 1.311 0.1951
M4 7.948 240.258 0.033 0.9737
M5 -144.623 240.479 -0.601 0.5499
M6 -326.973 243.015 -1.345 0.1837
M7 -181.870 239.794 -0.758 0.4513
M8 57.044 240.371 0.237 0.8132
M9 19.430 240.791 0.081 0.9360
M10 558.522 242.110 2.307 0.0247 *
M11 404.810 241.849 1.674 0.0996 .
t 37.845 5.620 6.734 8.18e-09 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 414.8 on 58 degrees of freedom
Multiple R-squared: 0.9439, Adjusted R-squared: 0.9313
F-statistic: 75.01 on 13 and 58 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.293393461 0.586786921 0.70660654
[2,] 0.161653887 0.323307775 0.83834611
[3,] 0.086586044 0.173172088 0.91341396
[4,] 0.048253666 0.096507331 0.95174633
[5,] 0.028226491 0.056452983 0.97177351
[6,] 0.012888812 0.025777623 0.98711119
[7,] 0.012138579 0.024277157 0.98786142
[8,] 0.017630274 0.035260548 0.98236973
[9,] 0.012308021 0.024616042 0.98769198
[10,] 0.007815083 0.015630166 0.99218492
[11,] 0.003611843 0.007223687 0.99638816
[12,] 0.001579372 0.003158743 0.99842063
[13,] 0.054711107 0.109422215 0.94528889
[14,] 0.056041039 0.112082079 0.94395896
[15,] 0.072388885 0.144777770 0.92761112
[16,] 0.176566699 0.353133399 0.82343330
[17,] 0.175204612 0.350409225 0.82479539
[18,] 0.178246394 0.356492788 0.82175361
[19,] 0.251609768 0.503219535 0.74839023
[20,] 0.275706168 0.551412337 0.72429383
[21,] 0.500676135 0.998647731 0.49932387
[22,] 0.463254300 0.926508599 0.53674570
[23,] 0.448740084 0.897480168 0.55125992
[24,] 0.466321415 0.932642830 0.53367859
[25,] 0.887606814 0.224786372 0.11239319
[26,] 0.910519468 0.178961063 0.08948053
[27,] 0.937104688 0.125790625 0.06289531
[28,] 0.944455537 0.111088926 0.05554446
[29,] 0.935265839 0.129468322 0.06473416
[30,] 0.967463959 0.065072083 0.03253604
[31,] 0.974938469 0.050123062 0.02506153
[32,] 0.972365811 0.055268379 0.02763419
[33,] 0.958381840 0.083236320 0.04161816
[34,] 0.929909326 0.140181348 0.07009067
[35,] 0.949123338 0.101753323 0.05087666
[36,] 0.910996195 0.178007609 0.08900380
[37,] 0.894441695 0.211116610 0.10555830
[38,] 0.834028396 0.331943207 0.16597160
[39,] 0.862205489 0.275589022 0.13779451
> postscript(file="/var/www/html/rcomp/tmp/15cdn1258564212.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/2nj4d1258564212.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/3pe7i1258564212.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/4609s1258564212.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/5tok51258564212.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
266.924968 209.974298 588.451032 870.491614 268.640893 302.471574
7 8 9 10 11 12
5.674377 -99.382803 61.684882 -41.204378 -313.136442 -334.918904
13 14 15 16 17 18
-53.591503 217.285102 25.162784 -147.446882 -344.122297 -349.364964
19 20 21 22 23 24
-468.688081 -749.395293 -494.656713 -520.943907 -597.199265 -198.485063
25 26 27 28 29 30
-612.378568 -153.954749 26.219983 -79.337842 245.542780 -17.206881
31 32 33 34 35 36
-138.877384 166.465307 -37.868170 -72.154934 134.161227 235.101500
37 38 39 40 41 42
-517.591746 -152.741599 -264.091280 -307.827210 -534.627167 247.852470
43 44 45 46 47 48
417.006231 641.421839 613.538567 398.028100 647.943314 687.959303
49 50 51 52 53 54
769.739706 3.714180 375.137783 -124.944241 129.903467 281.732642
55 56 57 58 59 60
120.613872 412.754367 224.246445 597.884826 284.825207 -119.057922
61 62 63 64 65 66
146.897142 -124.277232 -750.880302 -210.935438 234.662324 -465.484841
67 68 69 70 71 72
64.270985 -371.863417 -366.945010 -361.609707 -156.594041 -270.598914
> postscript(file="/var/www/html/rcomp/tmp/68pdn1258564212.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 266.924968 NA
1 209.974298 266.924968
2 588.451032 209.974298
3 870.491614 588.451032
4 268.640893 870.491614
5 302.471574 268.640893
6 5.674377 302.471574
7 -99.382803 5.674377
8 61.684882 -99.382803
9 -41.204378 61.684882
10 -313.136442 -41.204378
11 -334.918904 -313.136442
12 -53.591503 -334.918904
13 217.285102 -53.591503
14 25.162784 217.285102
15 -147.446882 25.162784
16 -344.122297 -147.446882
17 -349.364964 -344.122297
18 -468.688081 -349.364964
19 -749.395293 -468.688081
20 -494.656713 -749.395293
21 -520.943907 -494.656713
22 -597.199265 -520.943907
23 -198.485063 -597.199265
24 -612.378568 -198.485063
25 -153.954749 -612.378568
26 26.219983 -153.954749
27 -79.337842 26.219983
28 245.542780 -79.337842
29 -17.206881 245.542780
30 -138.877384 -17.206881
31 166.465307 -138.877384
32 -37.868170 166.465307
33 -72.154934 -37.868170
34 134.161227 -72.154934
35 235.101500 134.161227
36 -517.591746 235.101500
37 -152.741599 -517.591746
38 -264.091280 -152.741599
39 -307.827210 -264.091280
40 -534.627167 -307.827210
41 247.852470 -534.627167
42 417.006231 247.852470
43 641.421839 417.006231
44 613.538567 641.421839
45 398.028100 613.538567
46 647.943314 398.028100
47 687.959303 647.943314
48 769.739706 687.959303
49 3.714180 769.739706
50 375.137783 3.714180
51 -124.944241 375.137783
52 129.903467 -124.944241
53 281.732642 129.903467
54 120.613872 281.732642
55 412.754367 120.613872
56 224.246445 412.754367
57 597.884826 224.246445
58 284.825207 597.884826
59 -119.057922 284.825207
60 146.897142 -119.057922
61 -124.277232 146.897142
62 -750.880302 -124.277232
63 -210.935438 -750.880302
64 234.662324 -210.935438
65 -465.484841 234.662324
66 64.270985 -465.484841
67 -371.863417 64.270985
68 -366.945010 -371.863417
69 -361.609707 -366.945010
70 -156.594041 -361.609707
71 -270.598914 -156.594041
72 NA -270.598914
> dum1 <- dum[2:length(myerror),]
> dum1
lag(myerror, k = 1) myerror
[1,] 209.974298 266.924968
[2,] 588.451032 209.974298
[3,] 870.491614 588.451032
[4,] 268.640893 870.491614
[5,] 302.471574 268.640893
[6,] 5.674377 302.471574
[7,] -99.382803 5.674377
[8,] 61.684882 -99.382803
[9,] -41.204378 61.684882
[10,] -313.136442 -41.204378
[11,] -334.918904 -313.136442
[12,] -53.591503 -334.918904
[13,] 217.285102 -53.591503
[14,] 25.162784 217.285102
[15,] -147.446882 25.162784
[16,] -344.122297 -147.446882
[17,] -349.364964 -344.122297
[18,] -468.688081 -349.364964
[19,] -749.395293 -468.688081
[20,] -494.656713 -749.395293
[21,] -520.943907 -494.656713
[22,] -597.199265 -520.943907
[23,] -198.485063 -597.199265
[24,] -612.378568 -198.485063
[25,] -153.954749 -612.378568
[26,] 26.219983 -153.954749
[27,] -79.337842 26.219983
[28,] 245.542780 -79.337842
[29,] -17.206881 245.542780
[30,] -138.877384 -17.206881
[31,] 166.465307 -138.877384
[32,] -37.868170 166.465307
[33,] -72.154934 -37.868170
[34,] 134.161227 -72.154934
[35,] 235.101500 134.161227
[36,] -517.591746 235.101500
[37,] -152.741599 -517.591746
[38,] -264.091280 -152.741599
[39,] -307.827210 -264.091280
[40,] -534.627167 -307.827210
[41,] 247.852470 -534.627167
[42,] 417.006231 247.852470
[43,] 641.421839 417.006231
[44,] 613.538567 641.421839
[45,] 398.028100 613.538567
[46,] 647.943314 398.028100
[47,] 687.959303 647.943314
[48,] 769.739706 687.959303
[49,] 3.714180 769.739706
[50,] 375.137783 3.714180
[51,] -124.944241 375.137783
[52,] 129.903467 -124.944241
[53,] 281.732642 129.903467
[54,] 120.613872 281.732642
[55,] 412.754367 120.613872
[56,] 224.246445 412.754367
[57,] 597.884826 224.246445
[58,] 284.825207 597.884826
[59,] -119.057922 284.825207
[60,] 146.897142 -119.057922
[61,] -124.277232 146.897142
[62,] -750.880302 -124.277232
[63,] -210.935438 -750.880302
[64,] 234.662324 -210.935438
[65,] -465.484841 234.662324
[66,] 64.270985 -465.484841
[67,] -371.863417 64.270985
[68,] -366.945010 -371.863417
[69,] -361.609707 -366.945010
[70,] -156.594041 -361.609707
[71,] -270.598914 -156.594041
> z <- as.data.frame(dum1)
> z
lag(myerror, k = 1) myerror
1 209.974298 266.924968
2 588.451032 209.974298
3 870.491614 588.451032
4 268.640893 870.491614
5 302.471574 268.640893
6 5.674377 302.471574
7 -99.382803 5.674377
8 61.684882 -99.382803
9 -41.204378 61.684882
10 -313.136442 -41.204378
11 -334.918904 -313.136442
12 -53.591503 -334.918904
13 217.285102 -53.591503
14 25.162784 217.285102
15 -147.446882 25.162784
16 -344.122297 -147.446882
17 -349.364964 -344.122297
18 -468.688081 -349.364964
19 -749.395293 -468.688081
20 -494.656713 -749.395293
21 -520.943907 -494.656713
22 -597.199265 -520.943907
23 -198.485063 -597.199265
24 -612.378568 -198.485063
25 -153.954749 -612.378568
26 26.219983 -153.954749
27 -79.337842 26.219983
28 245.542780 -79.337842
29 -17.206881 245.542780
30 -138.877384 -17.206881
31 166.465307 -138.877384
32 -37.868170 166.465307
33 -72.154934 -37.868170
34 134.161227 -72.154934
35 235.101500 134.161227
36 -517.591746 235.101500
37 -152.741599 -517.591746
38 -264.091280 -152.741599
39 -307.827210 -264.091280
40 -534.627167 -307.827210
41 247.852470 -534.627167
42 417.006231 247.852470
43 641.421839 417.006231
44 613.538567 641.421839
45 398.028100 613.538567
46 647.943314 398.028100
47 687.959303 647.943314
48 769.739706 687.959303
49 3.714180 769.739706
50 375.137783 3.714180
51 -124.944241 375.137783
52 129.903467 -124.944241
53 281.732642 129.903467
54 120.613872 281.732642
55 412.754367 120.613872
56 224.246445 412.754367
57 597.884826 224.246445
58 284.825207 597.884826
59 -119.057922 284.825207
60 146.897142 -119.057922
61 -124.277232 146.897142
62 -750.880302 -124.277232
63 -210.935438 -750.880302
64 234.662324 -210.935438
65 -465.484841 234.662324
66 64.270985 -465.484841
67 -371.863417 64.270985
68 -366.945010 -371.863417
69 -361.609707 -366.945010
70 -156.594041 -361.609707
71 -270.598914 -156.594041
> 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/73x0u1258564212.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/8un4x1258564212.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/98clo1258564212.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/10v57f1258564212.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/11tqbw1258564212.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/1215vi1258564212.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/13qr7j1258564212.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/144dcs1258564212.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/15gi331258564213.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/16be7l1258564213.tab")
+ }
>
> system("convert tmp/15cdn1258564212.ps tmp/15cdn1258564212.png")
> system("convert tmp/2nj4d1258564212.ps tmp/2nj4d1258564212.png")
> system("convert tmp/3pe7i1258564212.ps tmp/3pe7i1258564212.png")
> system("convert tmp/4609s1258564212.ps tmp/4609s1258564212.png")
> system("convert tmp/5tok51258564212.ps tmp/5tok51258564212.png")
> system("convert tmp/68pdn1258564212.ps tmp/68pdn1258564212.png")
> system("convert tmp/73x0u1258564212.ps tmp/73x0u1258564212.png")
> system("convert tmp/8un4x1258564212.ps tmp/8un4x1258564212.png")
> system("convert tmp/98clo1258564212.ps tmp/98clo1258564212.png")
> system("convert tmp/10v57f1258564212.ps tmp/10v57f1258564212.png")
>
>
> proc.time()
user system elapsed
2.541 1.550 3.131