R version 2.15.2 (2012-10-26) -- "Trick or Treat"
Copyright (C) 2012 The R Foundation for Statistical Computing
ISBN 3-900051-07-0
Platform: i686-pc-linux-gnu (32-bit)
R is free software and comes with ABSOLUTELY NO WARRANTY.
You are welcome to redistribute it under certain conditions.
Type 'license()' or 'licence()' for distribution details.
R is a collaborative project with many contributors.
Type 'contributors()' for more information and
'citation()' on how to cite R or R packages in publications.
Type 'demo()' for some demos, 'help()' for on-line help, or
'help.start()' for an HTML browser interface to help.
Type 'q()' to quit R.
> x <- array(list(18308
+ ,4.041
+ ,79.6
+ ,7.2
+ ,1
+ ,1148
+ ,0.55
+ ,1
+ ,8.5
+ ,1
+ ,18068
+ ,3.665
+ ,32.3
+ ,5.7
+ ,1
+ ,7729
+ ,2.351
+ ,45.1
+ ,7.3
+ ,1
+ ,100484
+ ,29.76
+ ,190.8
+ ,7.5
+ ,1
+ ,16728
+ ,3.294
+ ,31.8
+ ,5
+ ,1
+ ,14630
+ ,3.287
+ ,678.4
+ ,6.7
+ ,1
+ ,4008
+ ,0.666
+ ,340.8
+ ,6.2
+ ,1
+ ,38927
+ ,12.938
+ ,239.6
+ ,7.3
+ ,1
+ ,22322
+ ,6.478
+ ,111.9
+ ,5
+ ,1
+ ,3711
+ ,1.108
+ ,172.5
+ ,2.8
+ ,1
+ ,3136
+ ,1.007
+ ,12.2
+ ,6.1
+ ,1
+ ,50508
+ ,11.431
+ ,205.6
+ ,7.1
+ ,1
+ ,28886
+ ,5.544
+ ,154.6
+ ,5.9
+ ,1
+ ,16996
+ ,2.777
+ ,49.7
+ ,4.6
+ ,1
+ ,13035
+ ,2.478
+ ,30.3
+ ,4.4
+ ,1
+ ,12973
+ ,3.685
+ ,92.8
+ ,7.4
+ ,1
+ ,16309
+ ,4.22
+ ,96.9
+ ,7.1
+ ,1
+ ,5227
+ ,1.228
+ ,39.8
+ ,7.5
+ ,1
+ ,19235
+ ,4.781
+ ,489.2
+ ,5.9
+ ,1
+ ,44487
+ ,6.016
+ ,767.6
+ ,9
+ ,1
+ ,44213
+ ,9.295
+ ,163.6
+ ,9.2
+ ,1
+ ,23619
+ ,4.375
+ ,55
+ ,5.1
+ ,1
+ ,9106
+ ,2.573
+ ,54.9
+ ,8.6
+ ,1
+ ,24917
+ ,5.117
+ ,74.3
+ ,6.6
+ ,1
+ ,3872
+ ,0.799
+ ,5.5
+ ,6.9
+ ,1
+ ,8945
+ ,1.578
+ ,20.5
+ ,2.7
+ ,1
+ ,2373
+ ,1.202
+ ,10.9
+ ,5.5
+ ,1
+ ,7128
+ ,1.109
+ ,123.7
+ ,7.2
+ ,1
+ ,23624
+ ,7.73
+ ,1042
+ ,6.6
+ ,1
+ ,5242
+ ,1.515
+ ,12.5
+ ,6.9
+ ,1
+ ,92629
+ ,17.99
+ ,381
+ ,7.2
+ ,1
+ ,28795
+ ,6.629
+ ,136.1
+ ,5.8
+ ,1
+ ,4487
+ ,0.639
+ ,9.3
+ ,4.1
+ ,1
+ ,48799
+ ,10.847
+ ,264.9
+ ,6.4
+ ,1
+ ,14067
+ ,3.146
+ ,45.8
+ ,6.7
+ ,1
+ ,12693
+ ,2.842
+ ,29.6
+ ,6
+ ,1
+ ,62184
+ ,11.882
+ ,265.1
+ ,6.9
+ ,1
+ ,9153
+ ,1.003
+ ,960.3
+ ,8.5
+ ,1
+ ,14250
+ ,3.487
+ ,115.8
+ ,6.2
+ ,1
+ ,3680
+ ,0.696
+ ,9.2
+ ,3.4
+ ,1
+ ,18063
+ ,4.877
+ ,118.3
+ ,6.6
+ ,1
+ ,65112
+ ,16.987
+ ,64.9
+ ,6.6
+ ,1
+ ,11340
+ ,1.723
+ ,21
+ ,4.9
+ ,1
+ ,4553
+ ,0.563
+ ,60.8
+ ,6.4
+ ,1
+ ,28960
+ ,6.187
+ ,156.3
+ ,5.8
+ ,1
+ ,19201
+ ,4.867
+ ,73.1
+ ,6.3
+ ,1
+ ,7533
+ ,1.793
+ ,74.5
+ ,10.5
+ ,1
+ ,26343
+ ,4.892
+ ,90.1
+ ,5.4
+ ,1
+ ,1641
+ ,0.454
+ ,4.7
+ ,5.1
+ ,1
+ ,145360
+ ,10.379
+ ,889
+ ,6.8
+ ,0
+ ,9066420
+ ,82.422
+ ,609
+ ,5.6
+ ,0
+ ,1038933
+ ,16.491
+ ,1259
+ ,3.8
+ ,0
+ ,2739420
+ ,60.876
+ ,289
+ ,8.2
+ ,0
+ ,61620
+ ,0.474
+ ,475
+ ,4.1
+ ,0
+ ,827530
+ ,7.523
+ ,490
+ ,2.8
+ ,0
+ ,534100
+ ,5.45
+ ,333
+ ,6.3
+ ,0
+ ,328755
+ ,10.605
+ ,300
+ ,11.4
+ ,0
+ ,1413895
+ ,40.397
+ ,210
+ ,19.4
+ ,0
+ ,2909136
+ ,60.607
+ ,650
+ ,5.8
+ ,0
+ ,3604246
+ ,58.133
+ ,512
+ ,6.9
+ ,0
+ ,917504
+ ,8.192
+ ,256
+ ,3.5
+ ,0)
+ ,dim=c(5
+ ,62)
+ ,dimnames=list(c('Bachelor'
+ ,'Populatie'
+ ,'Bevolkingsdicht'
+ ,'Werkloosheid'
+ ,'Land')
+ ,1:62))
> y <- array(NA,dim=c(5,62),dimnames=list(c('Bachelor','Populatie','Bevolkingsdicht','Werkloosheid','Land'),1:62))
> 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 = 'Do not include Seasonal Dummies'
> par1 = '1'
> library(lattice)
> library(lmtest)
Loading required package: zoo
Attaching package: 'zoo'
The following object(s) are masked from 'package:base':
as.Date, 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
Bachelor Populatie Bevolkingsdicht Werkloosheid Land
1 18308 4.041 79.6 7.2 1
2 1148 0.550 1.0 8.5 1
3 18068 3.665 32.3 5.7 1
4 7729 2.351 45.1 7.3 1
5 100484 29.760 190.8 7.5 1
6 16728 3.294 31.8 5.0 1
7 14630 3.287 678.4 6.7 1
8 4008 0.666 340.8 6.2 1
9 38927 12.938 239.6 7.3 1
10 22322 6.478 111.9 5.0 1
11 3711 1.108 172.5 2.8 1
12 3136 1.007 12.2 6.1 1
13 50508 11.431 205.6 7.1 1
14 28886 5.544 154.6 5.9 1
15 16996 2.777 49.7 4.6 1
16 13035 2.478 30.3 4.4 1
17 12973 3.685 92.8 7.4 1
18 16309 4.220 96.9 7.1 1
19 5227 1.228 39.8 7.5 1
20 19235 4.781 489.2 5.9 1
21 44487 6.016 767.6 9.0 1
22 44213 9.295 163.6 9.2 1
23 23619 4.375 55.0 5.1 1
24 9106 2.573 54.9 8.6 1
25 24917 5.117 74.3 6.6 1
26 3872 0.799 5.5 6.9 1
27 8945 1.578 20.5 2.7 1
28 2373 1.202 10.9 5.5 1
29 7128 1.109 123.7 7.2 1
30 23624 7.730 1042.0 6.6 1
31 5242 1.515 12.5 6.9 1
32 92629 17.990 381.0 7.2 1
33 28795 6.629 136.1 5.8 1
34 4487 0.639 9.3 4.1 1
35 48799 10.847 264.9 6.4 1
36 14067 3.146 45.8 6.7 1
37 12693 2.842 29.6 6.0 1
38 62184 11.882 265.1 6.9 1
39 9153 1.003 960.3 8.5 1
40 14250 3.487 115.8 6.2 1
41 3680 0.696 9.2 3.4 1
42 18063 4.877 118.3 6.6 1
43 65112 16.987 64.9 6.6 1
44 11340 1.723 21.0 4.9 1
45 4553 0.563 60.8 6.4 1
46 28960 6.187 156.3 5.8 1
47 19201 4.867 73.1 6.3 1
48 7533 1.793 74.5 10.5 1
49 26343 4.892 90.1 5.4 1
50 1641 0.454 4.7 5.1 1
51 145360 10.379 889.0 6.8 0
52 9066420 82.422 609.0 5.6 0
53 1038933 16.491 1259.0 3.8 0
54 2739420 60.876 289.0 8.2 0
55 61620 0.474 475.0 4.1 0
56 827530 7.523 490.0 2.8 0
57 534100 5.450 333.0 6.3 0
58 328755 10.605 300.0 11.4 0
59 1413895 40.397 210.0 19.4 0
60 2909136 60.607 650.0 5.8 0
61 3604246 58.133 512.0 6.9 0
62 917504 8.192 256.0 3.5 0
> k <- length(x[1,])
> df <- as.data.frame(x)
> (mylm <- lm(df))
Call:
lm(formula = df)
Coefficients:
(Intercept) Populatie Bevolkingsdicht Werkloosheid
574581.70 70010.42 -94.76 -94869.40
Land
-280628.86
> (mysum <- summary(mylm))
Call:
lm(formula = df)
Residuals:
Min 1Q Median 3Q Max
-1547379 -131954 8621 185810 3310414
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 574581.70 360351.53 1.595 0.1164
Populatie 70010.42 6087.17 11.501 <2e-16 ***
Bevolkingsdicht -94.76 315.30 -0.301 0.7649
Werkloosheid -94869.40 33513.86 -2.831 0.0064 **
Land -280628.86 267345.48 -1.050 0.2983
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 606800 on 57 degrees of freedom
Multiple R-squared: 0.8028, Adjusted R-squared: 0.7889
F-statistic: 58 on 4 and 57 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,] 9.168241e-08 1.833648e-07 0.9999999
[2,] 2.179673e-09 4.359346e-09 1.0000000
[3,] 2.293952e-11 4.587904e-11 1.0000000
[4,] 2.885348e-13 5.770695e-13 1.0000000
[5,] 1.893273e-15 3.786545e-15 1.0000000
[6,] 2.460180e-16 4.920360e-16 1.0000000
[7,] 5.561507e-18 1.112301e-17 1.0000000
[8,] 5.920559e-20 1.184112e-19 1.0000000
[9,] 4.515397e-22 9.030794e-22 1.0000000
[10,] 3.498171e-24 6.996342e-24 1.0000000
[11,] 2.407362e-26 4.814724e-26 1.0000000
[12,] 1.686813e-28 3.373625e-28 1.0000000
[13,] 1.056677e-30 2.113353e-30 1.0000000
[14,] 1.289077e-30 2.578155e-30 1.0000000
[15,] 2.505016e-32 5.010033e-32 1.0000000
[16,] 3.948201e-34 7.896402e-34 1.0000000
[17,] 4.226795e-36 8.453590e-36 1.0000000
[18,] 4.384729e-38 8.769458e-38 1.0000000
[19,] 3.745791e-40 7.491581e-40 1.0000000
[20,] 3.032654e-42 6.065308e-42 1.0000000
[21,] 2.835129e-44 5.670258e-44 1.0000000
[22,] 2.204698e-46 4.409396e-46 1.0000000
[23,] 2.075628e-47 4.151256e-47 1.0000000
[24,] 1.897904e-49 3.795809e-49 1.0000000
[25,] 2.328765e-48 4.657530e-48 1.0000000
[26,] 2.197390e-50 4.394780e-50 1.0000000
[27,] 1.890360e-52 3.780720e-52 1.0000000
[28,] 2.500505e-54 5.001010e-54 1.0000000
[29,] 2.040179e-56 4.080357e-56 1.0000000
[30,] 1.583601e-58 3.167202e-58 1.0000000
[31,] 1.117450e-59 2.234899e-59 1.0000000
[32,] 1.139867e-61 2.279735e-61 1.0000000
[33,] 8.594221e-64 1.718844e-63 1.0000000
[34,] 6.045030e-66 1.209006e-65 1.0000000
[35,] 4.463596e-68 8.927193e-68 1.0000000
[36,] 9.886041e-70 1.977208e-69 1.0000000
[37,] 6.978653e-72 1.395731e-71 1.0000000
[38,] 4.356311e-74 8.712621e-74 1.0000000
[39,] 3.097068e-76 6.194136e-76 1.0000000
[40,] 1.850254e-78 3.700507e-78 1.0000000
[41,] 1.789710e-80 3.579421e-80 1.0000000
[42,] 1.280929e-82 2.561859e-82 1.0000000
[43,] 5.885702e-85 1.177140e-84 1.0000000
[44,] 2.475967e-87 4.951934e-87 1.0000000
[45,] 7.736573e-01 4.526855e-01 0.2263427
[46,] 6.574423e-01 6.851154e-01 0.3425577
[47,] 8.459928e-01 3.080143e-01 0.1540072
> postscript(file="/var/wessaorg/rcomp/tmp/1m8k81356098369.ps",horizontal=F,onefile=F,pagecentre=F,paper="special",width=8.3333333333333,height=5.5555555555556)
> plot(x[,1], type='l', main='Actuals and Interpolation', ylab='value of Actuals and Interpolation (dots)', xlab='time or index')
> points(x[,1]-mysum$resid)
> grid()
> dev.off()
null device
1
> postscript(file="/var/wessaorg/rcomp/tmp/2v5gy1356098369.ps",horizontal=F,onefile=F,pagecentre=F,paper="special",width=8.3333333333333,height=5.5555555555556)
> plot(mysum$resid, type='b', pch=19, main='Residuals', ylab='value of Residuals', xlab='time or index')
> grid()
> dev.off()
null device
1
> postscript(file="/var/wessaorg/rcomp/tmp/3xo491356098369.ps",horizontal=F,onefile=F,pagecentre=F,paper="special",width=8.3333333333333,height=5.5555555555556)
> hist(mysum$resid, main='Residual Histogram', xlab='values of Residuals')
> grid()
> dev.off()
null device
1
> postscript(file="/var/wessaorg/rcomp/tmp/4auua1356098369.ps",horizontal=F,onefile=F,pagecentre=F,paper="special",width=8.3333333333333,height=5.5555555555556)
> densityplot(~mysum$resid,col='black',main='Residual Density Plot', xlab='values of Residuals')
> dev.off()
null device
1
> postscript(file="/var/wessaorg/rcomp/tmp/562ee1356098369.ps",horizontal=F,onefile=F,pagecentre=F,paper="special",width=8.3333333333333,height=5.5555555555556)
> qqnorm(mysum$resid, main='Residual Normal Q-Q Plot')
> qqline(mysum$resid)
> grid()
> dev.off()
null device
1
> (myerror <- as.ts(mysum$resid))
Time Series:
Start = 1
End = 62
Frequency = 1
1 2 3 4 5 6
132045.295 475174.097 11343.166 246001.774 -1547379.077 -30478.927
7 8 9 10 11 12
190460.180 283911.261 -445570.568 -240208.175 -85833.698 218542.036
13 14 15 16 17 18
-350679.438 -78825.911 -30267.172 -34107.197 171858.661 109666.763
19 20 21 22 23 24
340593.153 -3353.693 255910.574 -12186.176 -87583.919 356095.289
25 26 27 28 29 30
5899.235 309100.861 -137394.406 147082.181 330314.578 -86635.880
31 32 33 34 35 36
261006.690 -741649.698 -166118.139 55643.277 -372291.918 139826.178
37 38 39 40 41 42
91791.722 -383914.053 542363.544 75333.824 -15572.373 20016.988
43 44 45 46 47 48
-785820.171 63609.143 284109.612 -133094.467 -10888.687 591239.492
49 50 51 52 53 54
-89268.561 160182.730 -426510.096 3310414.230 -210389.376 -1291802.620
55 56 57 58 59 60
-112173.125 38324.529 207192.388 121650.672 -128532.609 -1296733.562
61 62
-337137.728 125697.298
> postscript(file="/var/wessaorg/rcomp/tmp/6rp5l1356098369.ps",horizontal=F,onefile=F,pagecentre=F,paper="special",width=8.3333333333333,height=5.5555555555556)
> dum <- cbind(lag(myerror,k=1),myerror)
> dum
Time Series:
Start = 0
End = 62
Frequency = 1
lag(myerror, k = 1) myerror
0 132045.295 NA
1 475174.097 132045.295
2 11343.166 475174.097
3 246001.774 11343.166
4 -1547379.077 246001.774
5 -30478.927 -1547379.077
6 190460.180 -30478.927
7 283911.261 190460.180
8 -445570.568 283911.261
9 -240208.175 -445570.568
10 -85833.698 -240208.175
11 218542.036 -85833.698
12 -350679.438 218542.036
13 -78825.911 -350679.438
14 -30267.172 -78825.911
15 -34107.197 -30267.172
16 171858.661 -34107.197
17 109666.763 171858.661
18 340593.153 109666.763
19 -3353.693 340593.153
20 255910.574 -3353.693
21 -12186.176 255910.574
22 -87583.919 -12186.176
23 356095.289 -87583.919
24 5899.235 356095.289
25 309100.861 5899.235
26 -137394.406 309100.861
27 147082.181 -137394.406
28 330314.578 147082.181
29 -86635.880 330314.578
30 261006.690 -86635.880
31 -741649.698 261006.690
32 -166118.139 -741649.698
33 55643.277 -166118.139
34 -372291.918 55643.277
35 139826.178 -372291.918
36 91791.722 139826.178
37 -383914.053 91791.722
38 542363.544 -383914.053
39 75333.824 542363.544
40 -15572.373 75333.824
41 20016.988 -15572.373
42 -785820.171 20016.988
43 63609.143 -785820.171
44 284109.612 63609.143
45 -133094.467 284109.612
46 -10888.687 -133094.467
47 591239.492 -10888.687
48 -89268.561 591239.492
49 160182.730 -89268.561
50 -426510.096 160182.730
51 3310414.230 -426510.096
52 -210389.376 3310414.230
53 -1291802.620 -210389.376
54 -112173.125 -1291802.620
55 38324.529 -112173.125
56 207192.388 38324.529
57 121650.672 207192.388
58 -128532.609 121650.672
59 -1296733.562 -128532.609
60 -337137.728 -1296733.562
61 125697.298 -337137.728
62 NA 125697.298
> dum1 <- dum[2:length(myerror),]
> dum1
lag(myerror, k = 1) myerror
[1,] 475174.097 132045.295
[2,] 11343.166 475174.097
[3,] 246001.774 11343.166
[4,] -1547379.077 246001.774
[5,] -30478.927 -1547379.077
[6,] 190460.180 -30478.927
[7,] 283911.261 190460.180
[8,] -445570.568 283911.261
[9,] -240208.175 -445570.568
[10,] -85833.698 -240208.175
[11,] 218542.036 -85833.698
[12,] -350679.438 218542.036
[13,] -78825.911 -350679.438
[14,] -30267.172 -78825.911
[15,] -34107.197 -30267.172
[16,] 171858.661 -34107.197
[17,] 109666.763 171858.661
[18,] 340593.153 109666.763
[19,] -3353.693 340593.153
[20,] 255910.574 -3353.693
[21,] -12186.176 255910.574
[22,] -87583.919 -12186.176
[23,] 356095.289 -87583.919
[24,] 5899.235 356095.289
[25,] 309100.861 5899.235
[26,] -137394.406 309100.861
[27,] 147082.181 -137394.406
[28,] 330314.578 147082.181
[29,] -86635.880 330314.578
[30,] 261006.690 -86635.880
[31,] -741649.698 261006.690
[32,] -166118.139 -741649.698
[33,] 55643.277 -166118.139
[34,] -372291.918 55643.277
[35,] 139826.178 -372291.918
[36,] 91791.722 139826.178
[37,] -383914.053 91791.722
[38,] 542363.544 -383914.053
[39,] 75333.824 542363.544
[40,] -15572.373 75333.824
[41,] 20016.988 -15572.373
[42,] -785820.171 20016.988
[43,] 63609.143 -785820.171
[44,] 284109.612 63609.143
[45,] -133094.467 284109.612
[46,] -10888.687 -133094.467
[47,] 591239.492 -10888.687
[48,] -89268.561 591239.492
[49,] 160182.730 -89268.561
[50,] -426510.096 160182.730
[51,] 3310414.230 -426510.096
[52,] -210389.376 3310414.230
[53,] -1291802.620 -210389.376
[54,] -112173.125 -1291802.620
[55,] 38324.529 -112173.125
[56,] 207192.388 38324.529
[57,] 121650.672 207192.388
[58,] -128532.609 121650.672
[59,] -1296733.562 -128532.609
[60,] -337137.728 -1296733.562
[61,] 125697.298 -337137.728
> z <- as.data.frame(dum1)
> z
lag(myerror, k = 1) myerror
1 475174.097 132045.295
2 11343.166 475174.097
3 246001.774 11343.166
4 -1547379.077 246001.774
5 -30478.927 -1547379.077
6 190460.180 -30478.927
7 283911.261 190460.180
8 -445570.568 283911.261
9 -240208.175 -445570.568
10 -85833.698 -240208.175
11 218542.036 -85833.698
12 -350679.438 218542.036
13 -78825.911 -350679.438
14 -30267.172 -78825.911
15 -34107.197 -30267.172
16 171858.661 -34107.197
17 109666.763 171858.661
18 340593.153 109666.763
19 -3353.693 340593.153
20 255910.574 -3353.693
21 -12186.176 255910.574
22 -87583.919 -12186.176
23 356095.289 -87583.919
24 5899.235 356095.289
25 309100.861 5899.235
26 -137394.406 309100.861
27 147082.181 -137394.406
28 330314.578 147082.181
29 -86635.880 330314.578
30 261006.690 -86635.880
31 -741649.698 261006.690
32 -166118.139 -741649.698
33 55643.277 -166118.139
34 -372291.918 55643.277
35 139826.178 -372291.918
36 91791.722 139826.178
37 -383914.053 91791.722
38 542363.544 -383914.053
39 75333.824 542363.544
40 -15572.373 75333.824
41 20016.988 -15572.373
42 -785820.171 20016.988
43 63609.143 -785820.171
44 284109.612 63609.143
45 -133094.467 284109.612
46 -10888.687 -133094.467
47 591239.492 -10888.687
48 -89268.561 591239.492
49 160182.730 -89268.561
50 -426510.096 160182.730
51 3310414.230 -426510.096
52 -210389.376 3310414.230
53 -1291802.620 -210389.376
54 -112173.125 -1291802.620
55 38324.529 -112173.125
56 207192.388 38324.529
57 121650.672 207192.388
58 -128532.609 121650.672
59 -1296733.562 -128532.609
60 -337137.728 -1296733.562
61 125697.298 -337137.728
> 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/wessaorg/rcomp/tmp/7w78b1356098369.ps",horizontal=F,onefile=F,pagecentre=F,paper="special",width=8.3333333333333,height=5.5555555555556)
> acf(mysum$resid, lag.max=length(mysum$resid)/2, main='Residual Autocorrelation Function')
> grid()
> dev.off()
null device
1
> postscript(file="/var/wessaorg/rcomp/tmp/8ue731356098369.ps",horizontal=F,onefile=F,pagecentre=F,paper="special",width=8.3333333333333,height=5.5555555555556)
> pacf(mysum$resid, lag.max=length(mysum$resid)/2, main='Residual Partial Autocorrelation Function')
> grid()
> dev.off()
null device
1
> postscript(file="/var/wessaorg/rcomp/tmp/9zg3g1356098369.ps",horizontal=F,onefile=F,pagecentre=F,paper="special",width=8.3333333333333,height=5.5555555555556)
> opar <- par(mfrow = c(2,2), oma = c(0, 0, 1.1, 0))
> plot(mylm, las = 1, sub='Residual Diagnostics')
> par(opar)
> dev.off()
null device
1
> if (n > n25) {
+ postscript(file="/var/wessaorg/rcomp/tmp/10s3wy1356098369.ps",horizontal=F,onefile=F,pagecentre=F,paper="special",width=8.3333333333333,height=5.5555555555556)
+ plot(kp3:nmkm3,gqarr[,2], main='Goldfeld-Quandt test',ylab='2-sided p-value',xlab='breakpoint')
+ grid()
+ dev.off()
+ }
null device
1
>
> #Note: the /var/wessaorg/rcomp/createtable file can be downloaded at http://www.wessa.net/cretab
> load(file="/var/wessaorg/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/wessaorg/rcomp/tmp/1105n21356098369.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/wessaorg/rcomp/tmp/12jljw1356098369.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/wessaorg/rcomp/tmp/13f32l1356098369.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/wessaorg/rcomp/tmp/14jzcy1356098369.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/wessaorg/rcomp/tmp/1525gm1356098369.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/wessaorg/rcomp/tmp/161i8g1356098369.tab")
+ }
>
> try(system("convert tmp/1m8k81356098369.ps tmp/1m8k81356098369.png",intern=TRUE))
character(0)
> try(system("convert tmp/2v5gy1356098369.ps tmp/2v5gy1356098369.png",intern=TRUE))
character(0)
> try(system("convert tmp/3xo491356098369.ps tmp/3xo491356098369.png",intern=TRUE))
character(0)
> try(system("convert tmp/4auua1356098369.ps tmp/4auua1356098369.png",intern=TRUE))
character(0)
> try(system("convert tmp/562ee1356098369.ps tmp/562ee1356098369.png",intern=TRUE))
character(0)
> try(system("convert tmp/6rp5l1356098369.ps tmp/6rp5l1356098369.png",intern=TRUE))
character(0)
> try(system("convert tmp/7w78b1356098369.ps tmp/7w78b1356098369.png",intern=TRUE))
character(0)
> try(system("convert tmp/8ue731356098369.ps tmp/8ue731356098369.png",intern=TRUE))
character(0)
> try(system("convert tmp/9zg3g1356098369.ps tmp/9zg3g1356098369.png",intern=TRUE))
character(0)
> try(system("convert tmp/10s3wy1356098369.ps tmp/10s3wy1356098369.png",intern=TRUE))
character(0)
>
>
> proc.time()
user system elapsed
6.044 0.906 7.356