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(277
+ ,52
+ ,99
+ ,104
+ ,172
+ ,79
+ ,8909
+ ,201
+ ,232
+ ,50
+ ,81
+ ,125
+ ,183
+ ,93
+ ,8841
+ ,201.2
+ ,256
+ ,59
+ ,95
+ ,98
+ ,162
+ ,68
+ ,8733
+ ,213.9
+ ,242
+ ,52
+ ,93
+ ,100
+ ,179
+ ,77
+ ,8885
+ ,209.7
+ ,302
+ ,66
+ ,109
+ ,93
+ ,162
+ ,78
+ ,8933
+ ,202.4
+ ,282
+ ,62
+ ,103
+ ,123
+ ,206
+ ,95
+ ,8854
+ ,187.8
+ ,288
+ ,59
+ ,101
+ ,116
+ ,194
+ ,88
+ ,8748
+ ,173.7
+ ,321
+ ,70
+ ,121
+ ,124
+ ,198
+ ,88
+ ,8827
+ ,172.3
+ ,316
+ ,74
+ ,112
+ ,126
+ ,219
+ ,102
+ ,8850
+ ,148
+ ,396
+ ,84
+ ,151
+ ,126
+ ,212
+ ,103
+ ,8761
+ ,129.8
+ ,362
+ ,71
+ ,142
+ ,156
+ ,265
+ ,131
+ ,8617
+ ,129.8
+ ,392
+ ,81
+ ,144
+ ,141
+ ,234
+ ,127
+ ,8758
+ ,117.9
+ ,414
+ ,92
+ ,154
+ ,163
+ ,259
+ ,133
+ ,8806
+ ,112.1
+ ,417
+ ,89
+ ,164
+ ,164
+ ,287
+ ,127
+ ,8710
+ ,94
+ ,476
+ ,100
+ ,188
+ ,156
+ ,278
+ ,138
+ ,8681
+ ,102.4
+ ,488
+ ,103
+ ,189
+ ,180
+ ,317
+ ,158
+ ,8819
+ ,115.8
+ ,489
+ ,97
+ ,188
+ ,187
+ ,320
+ ,167
+ ,8834
+ ,122.9
+ ,467
+ ,107
+ ,185
+ ,194
+ ,326
+ ,162
+ ,8742
+ ,120.9
+ ,460
+ ,93
+ ,188
+ ,168
+ ,316
+ ,149
+ ,8766
+ ,128.4
+ ,482
+ ,97
+ ,200
+ ,170
+ ,306
+ ,153
+ ,8902
+ ,148.8
+ ,510
+ ,100
+ ,211
+ ,177
+ ,315
+ ,166
+ ,8980
+ ,141.3
+ ,493
+ ,89
+ ,202
+ ,189
+ ,329
+ ,179
+ ,9031
+ ,163.7
+ ,476
+ ,102
+ ,198
+ ,194
+ ,316
+ ,176
+ ,8984
+ ,165.3
+ ,448
+ ,96
+ ,189
+ ,170
+ ,316
+ ,159
+ ,9150
+ ,187.3
+ ,410
+ ,81
+ ,174
+ ,156
+ ,297
+ ,151
+ ,9231
+ ,209.7
+ ,466
+ ,91
+ ,199
+ ,148
+ ,266
+ ,143
+ ,9230
+ ,230.1
+ ,417
+ ,84
+ ,175
+ ,167
+ ,296
+ ,169
+ ,9194
+ ,230.3
+ ,387
+ ,78
+ ,160
+ ,150
+ ,275
+ ,141
+ ,9307
+ ,234.9
+ ,370
+ ,70
+ ,160
+ ,141
+ ,252
+ ,134
+ ,9336
+ ,238.3
+ ,344
+ ,67
+ ,145
+ ,134
+ ,239
+ ,130
+ ,9310
+ ,246.8
+ ,396
+ ,76
+ ,172
+ ,127
+ ,231
+ ,112
+ ,9236
+ ,249.3
+ ,349
+ ,65
+ ,147
+ ,142
+ ,256
+ ,141
+ ,9244
+ ,247
+ ,326
+ ,66
+ ,138
+ ,132
+ ,232
+ ,116
+ ,9222
+ ,244.9
+ ,303
+ ,62
+ ,122
+ ,118
+ ,230
+ ,95
+ ,9186
+ ,246.7
+ ,300
+ ,66
+ ,118
+ ,115
+ ,205
+ ,98
+ ,9095
+ ,197.4
+ ,329
+ ,68
+ ,133
+ ,113
+ ,195
+ ,104
+ ,9216
+ ,153.9
+ ,304
+ ,59
+ ,118
+ ,123
+ ,207
+ ,121
+ ,9237
+ ,128.4
+ ,286
+ ,68
+ ,112
+ ,123
+ ,197
+ ,106
+ ,9207
+ ,130.7
+ ,281
+ ,68
+ ,109
+ ,103
+ ,194
+ ,90
+ ,9189
+ ,125.4
+ ,377
+ ,84
+ ,152
+ ,101
+ ,181
+ ,99
+ ,9183
+ ,115.6
+ ,344
+ ,75
+ ,141
+ ,135
+ ,246
+ ,130
+ ,9277
+ ,117.5
+ ,369
+ ,79
+ ,144
+ ,122
+ ,220
+ ,123
+ ,9305
+ ,125.3
+ ,390
+ ,92
+ ,152
+ ,142
+ ,234
+ ,133
+ ,9268
+ ,128.3
+ ,406
+ ,88
+ ,172
+ ,140
+ ,264
+ ,126
+ ,9259
+ ,134.7
+ ,426
+ ,98
+ ,168
+ ,138
+ ,266
+ ,137
+ ,9197
+ ,134.7
+ ,467
+ ,104
+ ,185
+ ,153
+ ,282
+ ,142
+ ,9293
+ ,134.1
+ ,437
+ ,95
+ ,174
+ ,172
+ ,312
+ ,153
+ ,9270
+ ,122.7
+ ,410
+ ,99
+ ,159
+ ,160
+ ,297
+ ,138
+ ,9395
+ ,117.8
+ ,390
+ ,93
+ ,155
+ ,146
+ ,269
+ ,139
+ ,9316
+ ,109.1
+ ,418
+ ,102
+ ,171
+ ,136
+ ,252
+ ,137
+ ,9237
+ ,108
+ ,398
+ ,91
+ ,161
+ ,139
+ ,265
+ ,152
+ ,9207
+ ,134.7
+ ,422
+ ,105
+ ,173
+ ,139
+ ,246
+ ,151
+ ,9189
+ ,134.7
+ ,439
+ ,100
+ ,179
+ ,140
+ ,263
+ ,158
+ ,9183
+ ,134.1
+ ,419
+ ,99
+ ,171
+ ,150
+ ,274
+ ,162
+ ,9277
+ ,122.7
+ ,484
+ ,111
+ ,202
+ ,142
+ ,262
+ ,156
+ ,9305
+ ,117.8
+ ,491
+ ,110
+ ,199
+ ,130
+ ,298
+ ,186
+ ,9268
+ ,109.1)
+ ,dim=c(8
+ ,56)
+ ,dimnames=list(c('werkeloosheid'
+ ,'onderwijshoog'
+ ,'onderwijsmiddelbaar'
+ ,'onderwijslaag'
+ ,'autochtoon'
+ ,'allochtonen'
+ ,'banen'
+ ,'vacatures')
+ ,1:56))
> y <- array(NA,dim=c(8,56),dimnames=list(c('werkeloosheid','onderwijshoog','onderwijsmiddelbaar','onderwijslaag','autochtoon','allochtonen','banen','vacatures'),1:56))
> 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'
> par3 <- 'No Linear Trend'
> par2 <- 'Do not include Seasonal 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, 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
werkeloosheid onderwijshoog onderwijsmiddelbaar onderwijslaag autochtoon
1 277 52 99 104 172
2 232 50 81 125 183
3 256 59 95 98 162
4 242 52 93 100 179
5 302 66 109 93 162
6 282 62 103 123 206
7 288 59 101 116 194
8 321 70 121 124 198
9 316 74 112 126 219
10 396 84 151 126 212
11 362 71 142 156 265
12 392 81 144 141 234
13 414 92 154 163 259
14 417 89 164 164 287
15 476 100 188 156 278
16 488 103 189 180 317
17 489 97 188 187 320
18 467 107 185 194 326
19 460 93 188 168 316
20 482 97 200 170 306
21 510 100 211 177 315
22 493 89 202 189 329
23 476 102 198 194 316
24 448 96 189 170 316
25 410 81 174 156 297
26 466 91 199 148 266
27 417 84 175 167 296
28 387 78 160 150 275
29 370 70 160 141 252
30 344 67 145 134 239
31 396 76 172 127 231
32 349 65 147 142 256
33 326 66 138 132 232
34 303 62 122 118 230
35 300 66 118 115 205
36 329 68 133 113 195
37 304 59 118 123 207
38 286 68 112 123 197
39 281 68 109 103 194
40 377 84 152 101 181
41 344 75 141 135 246
42 369 79 144 122 220
43 390 92 152 142 234
44 406 88 172 140 264
45 426 98 168 138 266
46 467 104 185 153 282
47 437 95 174 172 312
48 410 99 159 160 297
49 390 93 155 146 269
50 418 102 171 136 252
51 398 91 161 139 265
52 422 105 173 139 246
53 439 100 179 140 263
54 419 99 171 150 274
55 484 111 202 142 262
56 491 110 199 130 298
allochtonen banen vacatures
1 79 8909 201.0
2 93 8841 201.2
3 68 8733 213.9
4 77 8885 209.7
5 78 8933 202.4
6 95 8854 187.8
7 88 8748 173.7
8 88 8827 172.3
9 102 8850 148.0
10 103 8761 129.8
11 131 8617 129.8
12 127 8758 117.9
13 133 8806 112.1
14 127 8710 94.0
15 138 8681 102.4
16 158 8819 115.8
17 167 8834 122.9
18 162 8742 120.9
19 149 8766 128.4
20 153 8902 148.8
21 166 8980 141.3
22 179 9031 163.7
23 176 8984 165.3
24 159 9150 187.3
25 151 9231 209.7
26 143 9230 230.1
27 169 9194 230.3
28 141 9307 234.9
29 134 9336 238.3
30 130 9310 246.8
31 112 9236 249.3
32 141 9244 247.0
33 116 9222 244.9
34 95 9186 246.7
35 98 9095 197.4
36 104 9216 153.9
37 121 9237 128.4
38 106 9207 130.7
39 90 9189 125.4
40 99 9183 115.6
41 130 9277 117.5
42 123 9305 125.3
43 133 9268 128.3
44 126 9259 134.7
45 137 9197 134.7
46 142 9293 134.1
47 153 9270 122.7
48 138 9395 117.8
49 139 9316 109.1
50 137 9237 108.0
51 152 9207 134.7
52 151 9189 134.7
53 158 9183 134.1
54 162 9277 122.7
55 156 9305 117.8
56 186 9268 109.1
> k <- length(x[1,])
> df <- as.data.frame(x)
> (mylm <- lm(df))
Call:
lm(formula = df)
Coefficients:
(Intercept) onderwijshoog onderwijsmiddelbaar
336.94918 0.73647 1.76602
onderwijslaag autochtoon allochtonen
0.10237 -0.01110 0.04452
banen vacatures
-0.03193 -0.07557
> (mysum <- summary(mylm))
Call:
lm(formula = df)
Residuals:
Min 1Q Median 3Q Max
-12.427 -7.166 -1.030 6.489 16.943
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 336.949180 63.772687 5.284 3.05e-06 ***
onderwijshoog 0.736473 0.215405 3.419 0.00129 **
onderwijsmiddelbaar 1.766022 0.115506 15.289 < 2e-16 ***
onderwijslaag 0.102368 0.158366 0.646 0.52110
autochtoon -0.011102 0.104560 -0.106 0.91589
allochtonen 0.044522 0.125678 0.354 0.72470
banen -0.031931 0.007018 -4.550 3.66e-05 ***
vacatures -0.075575 0.038562 -1.960 0.05584 .
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 8.461 on 48 degrees of freedom
Multiple R-squared: 0.9886, Adjusted R-squared: 0.987
F-statistic: 597.1 on 7 and 48 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.6554810 0.6890381 0.34451905
[2,] 0.5385934 0.9228132 0.46140661
[3,] 0.4707905 0.9415809 0.52920954
[4,] 0.3491526 0.6983051 0.65084744
[5,] 0.2747319 0.5494638 0.72526811
[6,] 0.4052021 0.8104041 0.59479793
[7,] 0.5542953 0.8914094 0.44570470
[8,] 0.5196837 0.9606326 0.48031630
[9,] 0.4645888 0.9291777 0.53541117
[10,] 0.5537747 0.8924506 0.44622529
[11,] 0.5549972 0.8900056 0.44500280
[12,] 0.6035367 0.7929265 0.39646327
[13,] 0.7347924 0.5304151 0.26520756
[14,] 0.7483998 0.5032003 0.25160016
[15,] 0.7051792 0.5896416 0.29482079
[16,] 0.6461982 0.7076036 0.35380180
[17,] 0.5672738 0.8654524 0.43272618
[18,] 0.5049648 0.9900703 0.49503516
[19,] 0.5162329 0.9675341 0.48376705
[20,] 0.4487438 0.8974876 0.55125619
[21,] 0.3690010 0.7380019 0.63099905
[22,] 0.2876279 0.5752557 0.71237213
[23,] 0.2831354 0.5662708 0.71686458
[24,] 0.3127507 0.6255014 0.68724930
[25,] 0.2878518 0.5757035 0.71214825
[26,] 0.3368034 0.6736068 0.66319658
[27,] 0.3456013 0.6912026 0.65439869
[28,] 0.3224069 0.6448137 0.67759315
[29,] 0.2823261 0.5646522 0.71767392
[30,] 0.2741710 0.5483421 0.72582896
[31,] 0.2225186 0.4450372 0.77748139
[32,] 0.3303181 0.6606361 0.66968195
[33,] 0.7582760 0.4834481 0.24172403
[34,] 0.9696492 0.0607016 0.03035080
[35,] 0.9098812 0.1802375 0.09011877
> postscript(file="/var/wessaorg/rcomp/tmp/1lfv21353340584.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/2hm4l1353340584.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/3dc161353340584.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/4ffds1353340584.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/502tf1353340584.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 = 56
Frequency = 1
1 2 3 4 5 6
14.3288267 -2.2169471 -8.4144253 -9.6076631 13.2896932 -0.1336727
7 8 9 10 11 12
8.0525005 -0.7269238 5.5244280 4.9452812 -11.9137264 12.1619124
13 14 15 16 17 18
7.2531663 -9.1553171 -0.7029345 9.8265097 16.9429458 -10.6398829
19 20 21 22 23 24
-8.1648413 -4.9124832 2.1802274 10.8452037 -10.5672891 -8.0773471
25 26 27 28 29 30
-2.6822387 4.1432434 -1.2208374 6.3979251 -2.5497081 0.7125590
31 32 33 34 35 36
-2.3431001 0.4411628 -6.3920115 3.1426543 -2.4747042 -1.0352190
37 38 39 40 41 42
4.1794358 -10.0799840 -8.0308190 -1.0258300 -8.9653057 10.6280867
43 44 45 46 47 48
4.6338819 -10.6949378 6.7618924 14.7602196 7.1170656 9.0123015
49 50 51 52 53 54
-1.6070447 -10.1733454 -4.1826314 -12.4266807 -2.8026770 -6.8777022
55 56
-3.9854465 7.5025540
> postscript(file="/var/wessaorg/rcomp/tmp/6w4f61353340584.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 = 56
Frequency = 1
lag(myerror, k = 1) myerror
0 14.3288267 NA
1 -2.2169471 14.3288267
2 -8.4144253 -2.2169471
3 -9.6076631 -8.4144253
4 13.2896932 -9.6076631
5 -0.1336727 13.2896932
6 8.0525005 -0.1336727
7 -0.7269238 8.0525005
8 5.5244280 -0.7269238
9 4.9452812 5.5244280
10 -11.9137264 4.9452812
11 12.1619124 -11.9137264
12 7.2531663 12.1619124
13 -9.1553171 7.2531663
14 -0.7029345 -9.1553171
15 9.8265097 -0.7029345
16 16.9429458 9.8265097
17 -10.6398829 16.9429458
18 -8.1648413 -10.6398829
19 -4.9124832 -8.1648413
20 2.1802274 -4.9124832
21 10.8452037 2.1802274
22 -10.5672891 10.8452037
23 -8.0773471 -10.5672891
24 -2.6822387 -8.0773471
25 4.1432434 -2.6822387
26 -1.2208374 4.1432434
27 6.3979251 -1.2208374
28 -2.5497081 6.3979251
29 0.7125590 -2.5497081
30 -2.3431001 0.7125590
31 0.4411628 -2.3431001
32 -6.3920115 0.4411628
33 3.1426543 -6.3920115
34 -2.4747042 3.1426543
35 -1.0352190 -2.4747042
36 4.1794358 -1.0352190
37 -10.0799840 4.1794358
38 -8.0308190 -10.0799840
39 -1.0258300 -8.0308190
40 -8.9653057 -1.0258300
41 10.6280867 -8.9653057
42 4.6338819 10.6280867
43 -10.6949378 4.6338819
44 6.7618924 -10.6949378
45 14.7602196 6.7618924
46 7.1170656 14.7602196
47 9.0123015 7.1170656
48 -1.6070447 9.0123015
49 -10.1733454 -1.6070447
50 -4.1826314 -10.1733454
51 -12.4266807 -4.1826314
52 -2.8026770 -12.4266807
53 -6.8777022 -2.8026770
54 -3.9854465 -6.8777022
55 7.5025540 -3.9854465
56 NA 7.5025540
> dum1 <- dum[2:length(myerror),]
> dum1
lag(myerror, k = 1) myerror
[1,] -2.2169471 14.3288267
[2,] -8.4144253 -2.2169471
[3,] -9.6076631 -8.4144253
[4,] 13.2896932 -9.6076631
[5,] -0.1336727 13.2896932
[6,] 8.0525005 -0.1336727
[7,] -0.7269238 8.0525005
[8,] 5.5244280 -0.7269238
[9,] 4.9452812 5.5244280
[10,] -11.9137264 4.9452812
[11,] 12.1619124 -11.9137264
[12,] 7.2531663 12.1619124
[13,] -9.1553171 7.2531663
[14,] -0.7029345 -9.1553171
[15,] 9.8265097 -0.7029345
[16,] 16.9429458 9.8265097
[17,] -10.6398829 16.9429458
[18,] -8.1648413 -10.6398829
[19,] -4.9124832 -8.1648413
[20,] 2.1802274 -4.9124832
[21,] 10.8452037 2.1802274
[22,] -10.5672891 10.8452037
[23,] -8.0773471 -10.5672891
[24,] -2.6822387 -8.0773471
[25,] 4.1432434 -2.6822387
[26,] -1.2208374 4.1432434
[27,] 6.3979251 -1.2208374
[28,] -2.5497081 6.3979251
[29,] 0.7125590 -2.5497081
[30,] -2.3431001 0.7125590
[31,] 0.4411628 -2.3431001
[32,] -6.3920115 0.4411628
[33,] 3.1426543 -6.3920115
[34,] -2.4747042 3.1426543
[35,] -1.0352190 -2.4747042
[36,] 4.1794358 -1.0352190
[37,] -10.0799840 4.1794358
[38,] -8.0308190 -10.0799840
[39,] -1.0258300 -8.0308190
[40,] -8.9653057 -1.0258300
[41,] 10.6280867 -8.9653057
[42,] 4.6338819 10.6280867
[43,] -10.6949378 4.6338819
[44,] 6.7618924 -10.6949378
[45,] 14.7602196 6.7618924
[46,] 7.1170656 14.7602196
[47,] 9.0123015 7.1170656
[48,] -1.6070447 9.0123015
[49,] -10.1733454 -1.6070447
[50,] -4.1826314 -10.1733454
[51,] -12.4266807 -4.1826314
[52,] -2.8026770 -12.4266807
[53,] -6.8777022 -2.8026770
[54,] -3.9854465 -6.8777022
[55,] 7.5025540 -3.9854465
> z <- as.data.frame(dum1)
> z
lag(myerror, k = 1) myerror
1 -2.2169471 14.3288267
2 -8.4144253 -2.2169471
3 -9.6076631 -8.4144253
4 13.2896932 -9.6076631
5 -0.1336727 13.2896932
6 8.0525005 -0.1336727
7 -0.7269238 8.0525005
8 5.5244280 -0.7269238
9 4.9452812 5.5244280
10 -11.9137264 4.9452812
11 12.1619124 -11.9137264
12 7.2531663 12.1619124
13 -9.1553171 7.2531663
14 -0.7029345 -9.1553171
15 9.8265097 -0.7029345
16 16.9429458 9.8265097
17 -10.6398829 16.9429458
18 -8.1648413 -10.6398829
19 -4.9124832 -8.1648413
20 2.1802274 -4.9124832
21 10.8452037 2.1802274
22 -10.5672891 10.8452037
23 -8.0773471 -10.5672891
24 -2.6822387 -8.0773471
25 4.1432434 -2.6822387
26 -1.2208374 4.1432434
27 6.3979251 -1.2208374
28 -2.5497081 6.3979251
29 0.7125590 -2.5497081
30 -2.3431001 0.7125590
31 0.4411628 -2.3431001
32 -6.3920115 0.4411628
33 3.1426543 -6.3920115
34 -2.4747042 3.1426543
35 -1.0352190 -2.4747042
36 4.1794358 -1.0352190
37 -10.0799840 4.1794358
38 -8.0308190 -10.0799840
39 -1.0258300 -8.0308190
40 -8.9653057 -1.0258300
41 10.6280867 -8.9653057
42 4.6338819 10.6280867
43 -10.6949378 4.6338819
44 6.7618924 -10.6949378
45 14.7602196 6.7618924
46 7.1170656 14.7602196
47 9.0123015 7.1170656
48 -1.6070447 9.0123015
49 -10.1733454 -1.6070447
50 -4.1826314 -10.1733454
51 -12.4266807 -4.1826314
52 -2.8026770 -12.4266807
53 -6.8777022 -2.8026770
54 -3.9854465 -6.8777022
55 7.5025540 -3.9854465
> 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/7mow51353340584.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/8pv311353340584.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/9vyoq1353340584.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/10pxkg1353340584.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/112ra91353340584.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/126zmz1353340584.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/134oav1353340584.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/14dsmy1353340584.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/15p77v1353340584.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/16ab0l1353340584.tab")
+ }
>
> try(system("convert tmp/1lfv21353340584.ps tmp/1lfv21353340584.png",intern=TRUE))
character(0)
> try(system("convert tmp/2hm4l1353340584.ps tmp/2hm4l1353340584.png",intern=TRUE))
character(0)
> try(system("convert tmp/3dc161353340584.ps tmp/3dc161353340584.png",intern=TRUE))
character(0)
> try(system("convert tmp/4ffds1353340584.ps tmp/4ffds1353340584.png",intern=TRUE))
character(0)
> try(system("convert tmp/502tf1353340584.ps tmp/502tf1353340584.png",intern=TRUE))
character(0)
> try(system("convert tmp/6w4f61353340584.ps tmp/6w4f61353340584.png",intern=TRUE))
character(0)
> try(system("convert tmp/7mow51353340584.ps tmp/7mow51353340584.png",intern=TRUE))
character(0)
> try(system("convert tmp/8pv311353340584.ps tmp/8pv311353340584.png",intern=TRUE))
character(0)
> try(system("convert tmp/9vyoq1353340584.ps tmp/9vyoq1353340584.png",intern=TRUE))
character(0)
> try(system("convert tmp/10pxkg1353340584.ps tmp/10pxkg1353340584.png",intern=TRUE))
character(0)
>
>
> proc.time()
user system elapsed
6.119 0.840 7.170