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
+ ,232
+ ,50
+ ,81
+ ,125
+ ,183
+ ,93
+ ,8841
+ ,256
+ ,59
+ ,95
+ ,98
+ ,162
+ ,68
+ ,8733
+ ,242
+ ,52
+ ,93
+ ,100
+ ,179
+ ,77
+ ,8885
+ ,302
+ ,66
+ ,109
+ ,93
+ ,162
+ ,78
+ ,8933
+ ,282
+ ,62
+ ,103
+ ,123
+ ,206
+ ,95
+ ,8854
+ ,288
+ ,59
+ ,101
+ ,116
+ ,194
+ ,88
+ ,8748
+ ,321
+ ,70
+ ,121
+ ,124
+ ,198
+ ,88
+ ,8827
+ ,316
+ ,74
+ ,112
+ ,126
+ ,219
+ ,102
+ ,8850
+ ,396
+ ,84
+ ,151
+ ,126
+ ,212
+ ,103
+ ,8761
+ ,362
+ ,71
+ ,142
+ ,156
+ ,265
+ ,131
+ ,8617
+ ,392
+ ,81
+ ,144
+ ,141
+ ,234
+ ,127
+ ,8758
+ ,414
+ ,92
+ ,154
+ ,163
+ ,259
+ ,133
+ ,8806
+ ,417
+ ,89
+ ,164
+ ,164
+ ,287
+ ,127
+ ,8710
+ ,476
+ ,100
+ ,188
+ ,156
+ ,278
+ ,138
+ ,8681
+ ,488
+ ,103
+ ,189
+ ,180
+ ,317
+ ,158
+ ,8819
+ ,489
+ ,97
+ ,188
+ ,187
+ ,320
+ ,167
+ ,8834
+ ,467
+ ,107
+ ,185
+ ,194
+ ,326
+ ,162
+ ,8742
+ ,460
+ ,93
+ ,188
+ ,168
+ ,316
+ ,149
+ ,8766
+ ,482
+ ,97
+ ,200
+ ,170
+ ,306
+ ,153
+ ,8902
+ ,510
+ ,100
+ ,211
+ ,177
+ ,315
+ ,166
+ ,8980
+ ,493
+ ,89
+ ,202
+ ,189
+ ,329
+ ,179
+ ,9031
+ ,476
+ ,102
+ ,198
+ ,194
+ ,316
+ ,176
+ ,8984
+ ,448
+ ,96
+ ,189
+ ,170
+ ,316
+ ,159
+ ,9150
+ ,410
+ ,81
+ ,174
+ ,156
+ ,297
+ ,151
+ ,9231
+ ,466
+ ,91
+ ,199
+ ,148
+ ,266
+ ,143
+ ,9230
+ ,417
+ ,84
+ ,175
+ ,167
+ ,296
+ ,169
+ ,9194
+ ,387
+ ,78
+ ,160
+ ,150
+ ,275
+ ,141
+ ,9307
+ ,370
+ ,70
+ ,160
+ ,141
+ ,252
+ ,134
+ ,9336
+ ,344
+ ,67
+ ,145
+ ,134
+ ,239
+ ,130
+ ,9310
+ ,396
+ ,76
+ ,172
+ ,127
+ ,231
+ ,112
+ ,9236
+ ,349
+ ,65
+ ,147
+ ,142
+ ,256
+ ,141
+ ,9244
+ ,326
+ ,66
+ ,138
+ ,132
+ ,232
+ ,116
+ ,9222
+ ,303
+ ,62
+ ,122
+ ,118
+ ,230
+ ,95
+ ,9186
+ ,300
+ ,66
+ ,118
+ ,115
+ ,205
+ ,98
+ ,9095
+ ,329
+ ,68
+ ,133
+ ,113
+ ,195
+ ,104
+ ,9216
+ ,304
+ ,59
+ ,118
+ ,123
+ ,207
+ ,121
+ ,9237
+ ,286
+ ,68
+ ,112
+ ,123
+ ,197
+ ,106
+ ,9207
+ ,281
+ ,68
+ ,109
+ ,103
+ ,194
+ ,90
+ ,9189
+ ,377
+ ,84
+ ,152
+ ,101
+ ,181
+ ,99
+ ,9183
+ ,344
+ ,75
+ ,141
+ ,135
+ ,246
+ ,130
+ ,9277
+ ,369
+ ,79
+ ,144
+ ,122
+ ,220
+ ,123
+ ,9305
+ ,390
+ ,92
+ ,152
+ ,142
+ ,234
+ ,133
+ ,9268
+ ,406
+ ,88
+ ,172
+ ,140
+ ,264
+ ,126
+ ,9259
+ ,426
+ ,98
+ ,168
+ ,138
+ ,266
+ ,137
+ ,9197
+ ,467
+ ,104
+ ,185
+ ,153
+ ,282
+ ,142
+ ,9293
+ ,437
+ ,95
+ ,174
+ ,172
+ ,312
+ ,153
+ ,9270
+ ,410
+ ,99
+ ,159
+ ,160
+ ,297
+ ,138
+ ,9395
+ ,390
+ ,93
+ ,155
+ ,146
+ ,269
+ ,139
+ ,9316
+ ,418
+ ,102
+ ,171
+ ,136
+ ,252
+ ,137
+ ,9237
+ ,398
+ ,91
+ ,161
+ ,139
+ ,265
+ ,152
+ ,9207
+ ,422
+ ,105
+ ,173
+ ,139
+ ,246
+ ,151
+ ,9189
+ ,439
+ ,100
+ ,179
+ ,140
+ ,263
+ ,158
+ ,9183
+ ,419
+ ,99
+ ,171
+ ,150
+ ,274
+ ,162
+ ,9277
+ ,484
+ ,111
+ ,202
+ ,142
+ ,262
+ ,156
+ ,9305
+ ,491
+ ,110
+ ,199
+ ,130
+ ,298
+ ,186
+ ,9268)
+ ,dim=c(7
+ ,56)
+ ,dimnames=list(c('werkeloosheid'
+ ,'onderwijshoog'
+ ,'onderwijsmiddelbaar'
+ ,'onderwijslaag'
+ ,'autochtoon'
+ ,'allochtonen'
+ ,'banen')
+ ,1:56))
> y <- array(NA,dim=c(7,56),dimnames=list(c('werkeloosheid','onderwijshoog','onderwijsmiddelbaar','onderwijslaag','autochtoon','allochtonen','banen'),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
1 79 8909
2 93 8841
3 68 8733
4 77 8885
5 78 8933
6 95 8854
7 88 8748
8 88 8827
9 102 8850
10 103 8761
11 131 8617
12 127 8758
13 133 8806
14 127 8710
15 138 8681
16 158 8819
17 167 8834
18 162 8742
19 149 8766
20 153 8902
21 166 8980
22 179 9031
23 176 8984
24 159 9150
25 151 9231
26 143 9230
27 169 9194
28 141 9307
29 134 9336
30 130 9310
31 112 9236
32 141 9244
33 116 9222
34 95 9186
35 98 9095
36 104 9216
37 121 9237
38 106 9207
39 90 9189
40 99 9183
41 130 9277
42 123 9305
43 133 9268
44 126 9259
45 137 9197
46 142 9293
47 153 9270
48 138 9395
49 139 9316
50 137 9237
51 152 9207
52 151 9189
53 158 9183
54 162 9277
55 156 9305
56 186 9268
> k <- length(x[1,])
> df <- as.data.frame(x)
> (mylm <- lm(df))
Call:
lm(formula = df)
Coefficients:
(Intercept) onderwijshoog onderwijsmiddelbaar
347.57759 1.04313 1.65856
onderwijslaag autochtoon allochtonen
0.09636 -0.01959 0.07299
banen
-0.03548
> (mysum <- summary(mylm))
Call:
lm(formula = df)
Residuals:
Min 1Q Median 3Q Max
-15.633 -5.205 -1.500 6.094 17.896
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 347.577587 65.357759 5.318 2.58e-06 ***
onderwijshoog 1.043132 0.152271 6.850 1.13e-08 ***
onderwijsmiddelbaar 1.658560 0.104567 15.861 < 2e-16 ***
onderwijslaag 0.096357 0.162861 0.592 0.557
autochtoon -0.019589 0.107456 -0.182 0.856
allochtonen 0.072989 0.128404 0.568 0.572
banen -0.035477 0.006974 -5.087 5.74e-06 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 8.703 on 49 degrees of freedom
Multiple R-squared: 0.9877, Adjusted R-squared: 0.9862
F-statistic: 657.8 on 6 and 49 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.8155788 0.3688423 0.18442115
[2,] 0.7736456 0.4527088 0.22635440
[3,] 0.6984094 0.6031812 0.30159061
[4,] 0.6017597 0.7964805 0.39824026
[5,] 0.4843217 0.9686435 0.51567826
[6,] 0.4060222 0.8120444 0.59397781
[7,] 0.3372208 0.6744416 0.66277918
[8,] 0.4233041 0.8466082 0.57669592
[9,] 0.5600930 0.8798140 0.43990701
[10,] 0.5262310 0.9475380 0.47376898
[11,] 0.6550298 0.6899405 0.34497024
[12,] 0.6578763 0.6842474 0.34212371
[13,] 0.8044663 0.3910674 0.19553368
[14,] 0.9481176 0.1037647 0.05188237
[15,] 0.9476272 0.1047457 0.05237284
[16,] 0.9225522 0.1548956 0.07744778
[17,] 0.8979596 0.2040809 0.10204044
[18,] 0.8843872 0.2312255 0.11561277
[19,] 0.8511440 0.2977119 0.14885595
[20,] 0.8159391 0.3681218 0.18406088
[21,] 0.7694578 0.4610844 0.23054221
[22,] 0.7134862 0.5730275 0.28651375
[23,] 0.6415091 0.7169819 0.35849093
[24,] 0.6670694 0.6658612 0.33293062
[25,] 0.6308901 0.7382197 0.36910987
[26,] 0.5488320 0.9023359 0.45116795
[27,] 0.4534795 0.9069590 0.54652048
[28,] 0.4541144 0.9082287 0.54588563
[29,] 0.3789374 0.7578749 0.62106257
[30,] 0.3135440 0.6270880 0.68645598
[31,] 0.2526536 0.5053071 0.74734644
[32,] 0.2352824 0.4705648 0.76471758
[33,] 0.3675741 0.7351481 0.63242593
[34,] 0.7941925 0.4116149 0.20580746
[35,] 0.9079548 0.1840904 0.09204518
[36,] 0.8294967 0.3410066 0.17050330
[37,] 0.9440280 0.1119440 0.05597202
> postscript(file="/var/wessaorg/rcomp/tmp/1zumd1353002011.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/2s1cq1353002011.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/3fx0j1353002011.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/4378u1353002011.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/5sls31353002011.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.6313449 -3.6706182 -12.0952063 -10.6002199 10.2303722 -1.7180453
7 8 9 10 11 12
7.9182317 -1.6172175 4.1500813 5.6673136 -8.8498434 13.5340647
13 14 15 16 17 18
7.1088647 -5.8630952 0.6198309 9.7193668 17.8962441 -13.0153180
19 20 21 22 23 24
-6.2774584 -4.2083592 3.7382522 13.1181931 -12.9931811 -10.3647400
25 26 27 28 29 30
-3.4049794 1.4117235 -4.8989141 3.5176131 -3.1809029 -1.3837131
31 32 33 34 35 36
-4.3467592 -1.1967951 -8.7752182 0.4996840 -3.6866568 2.2003124
37 38 39 40 41 42
10.2426116 -7.3595789 -4.9862955 2.0737843 -4.2245097 13.8748939
43 44 45 46 47 48
4.3502587 -7.6764131 5.7559271 14.2105570 8.9809418 9.0787685
49 50 51 52 53 54
-0.1033938 -10.0547217 -4.1882402 -15.6326022 -3.8554403 -7.2490033
55 56
-4.2148742 7.1630781
> postscript(file="/var/wessaorg/rcomp/tmp/6065q1353002011.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.6313449 NA
1 -3.6706182 14.6313449
2 -12.0952063 -3.6706182
3 -10.6002199 -12.0952063
4 10.2303722 -10.6002199
5 -1.7180453 10.2303722
6 7.9182317 -1.7180453
7 -1.6172175 7.9182317
8 4.1500813 -1.6172175
9 5.6673136 4.1500813
10 -8.8498434 5.6673136
11 13.5340647 -8.8498434
12 7.1088647 13.5340647
13 -5.8630952 7.1088647
14 0.6198309 -5.8630952
15 9.7193668 0.6198309
16 17.8962441 9.7193668
17 -13.0153180 17.8962441
18 -6.2774584 -13.0153180
19 -4.2083592 -6.2774584
20 3.7382522 -4.2083592
21 13.1181931 3.7382522
22 -12.9931811 13.1181931
23 -10.3647400 -12.9931811
24 -3.4049794 -10.3647400
25 1.4117235 -3.4049794
26 -4.8989141 1.4117235
27 3.5176131 -4.8989141
28 -3.1809029 3.5176131
29 -1.3837131 -3.1809029
30 -4.3467592 -1.3837131
31 -1.1967951 -4.3467592
32 -8.7752182 -1.1967951
33 0.4996840 -8.7752182
34 -3.6866568 0.4996840
35 2.2003124 -3.6866568
36 10.2426116 2.2003124
37 -7.3595789 10.2426116
38 -4.9862955 -7.3595789
39 2.0737843 -4.9862955
40 -4.2245097 2.0737843
41 13.8748939 -4.2245097
42 4.3502587 13.8748939
43 -7.6764131 4.3502587
44 5.7559271 -7.6764131
45 14.2105570 5.7559271
46 8.9809418 14.2105570
47 9.0787685 8.9809418
48 -0.1033938 9.0787685
49 -10.0547217 -0.1033938
50 -4.1882402 -10.0547217
51 -15.6326022 -4.1882402
52 -3.8554403 -15.6326022
53 -7.2490033 -3.8554403
54 -4.2148742 -7.2490033
55 7.1630781 -4.2148742
56 NA 7.1630781
> dum1 <- dum[2:length(myerror),]
> dum1
lag(myerror, k = 1) myerror
[1,] -3.6706182 14.6313449
[2,] -12.0952063 -3.6706182
[3,] -10.6002199 -12.0952063
[4,] 10.2303722 -10.6002199
[5,] -1.7180453 10.2303722
[6,] 7.9182317 -1.7180453
[7,] -1.6172175 7.9182317
[8,] 4.1500813 -1.6172175
[9,] 5.6673136 4.1500813
[10,] -8.8498434 5.6673136
[11,] 13.5340647 -8.8498434
[12,] 7.1088647 13.5340647
[13,] -5.8630952 7.1088647
[14,] 0.6198309 -5.8630952
[15,] 9.7193668 0.6198309
[16,] 17.8962441 9.7193668
[17,] -13.0153180 17.8962441
[18,] -6.2774584 -13.0153180
[19,] -4.2083592 -6.2774584
[20,] 3.7382522 -4.2083592
[21,] 13.1181931 3.7382522
[22,] -12.9931811 13.1181931
[23,] -10.3647400 -12.9931811
[24,] -3.4049794 -10.3647400
[25,] 1.4117235 -3.4049794
[26,] -4.8989141 1.4117235
[27,] 3.5176131 -4.8989141
[28,] -3.1809029 3.5176131
[29,] -1.3837131 -3.1809029
[30,] -4.3467592 -1.3837131
[31,] -1.1967951 -4.3467592
[32,] -8.7752182 -1.1967951
[33,] 0.4996840 -8.7752182
[34,] -3.6866568 0.4996840
[35,] 2.2003124 -3.6866568
[36,] 10.2426116 2.2003124
[37,] -7.3595789 10.2426116
[38,] -4.9862955 -7.3595789
[39,] 2.0737843 -4.9862955
[40,] -4.2245097 2.0737843
[41,] 13.8748939 -4.2245097
[42,] 4.3502587 13.8748939
[43,] -7.6764131 4.3502587
[44,] 5.7559271 -7.6764131
[45,] 14.2105570 5.7559271
[46,] 8.9809418 14.2105570
[47,] 9.0787685 8.9809418
[48,] -0.1033938 9.0787685
[49,] -10.0547217 -0.1033938
[50,] -4.1882402 -10.0547217
[51,] -15.6326022 -4.1882402
[52,] -3.8554403 -15.6326022
[53,] -7.2490033 -3.8554403
[54,] -4.2148742 -7.2490033
[55,] 7.1630781 -4.2148742
> z <- as.data.frame(dum1)
> z
lag(myerror, k = 1) myerror
1 -3.6706182 14.6313449
2 -12.0952063 -3.6706182
3 -10.6002199 -12.0952063
4 10.2303722 -10.6002199
5 -1.7180453 10.2303722
6 7.9182317 -1.7180453
7 -1.6172175 7.9182317
8 4.1500813 -1.6172175
9 5.6673136 4.1500813
10 -8.8498434 5.6673136
11 13.5340647 -8.8498434
12 7.1088647 13.5340647
13 -5.8630952 7.1088647
14 0.6198309 -5.8630952
15 9.7193668 0.6198309
16 17.8962441 9.7193668
17 -13.0153180 17.8962441
18 -6.2774584 -13.0153180
19 -4.2083592 -6.2774584
20 3.7382522 -4.2083592
21 13.1181931 3.7382522
22 -12.9931811 13.1181931
23 -10.3647400 -12.9931811
24 -3.4049794 -10.3647400
25 1.4117235 -3.4049794
26 -4.8989141 1.4117235
27 3.5176131 -4.8989141
28 -3.1809029 3.5176131
29 -1.3837131 -3.1809029
30 -4.3467592 -1.3837131
31 -1.1967951 -4.3467592
32 -8.7752182 -1.1967951
33 0.4996840 -8.7752182
34 -3.6866568 0.4996840
35 2.2003124 -3.6866568
36 10.2426116 2.2003124
37 -7.3595789 10.2426116
38 -4.9862955 -7.3595789
39 2.0737843 -4.9862955
40 -4.2245097 2.0737843
41 13.8748939 -4.2245097
42 4.3502587 13.8748939
43 -7.6764131 4.3502587
44 5.7559271 -7.6764131
45 14.2105570 5.7559271
46 8.9809418 14.2105570
47 9.0787685 8.9809418
48 -0.1033938 9.0787685
49 -10.0547217 -0.1033938
50 -4.1882402 -10.0547217
51 -15.6326022 -4.1882402
52 -3.8554403 -15.6326022
53 -7.2490033 -3.8554403
54 -4.2148742 -7.2490033
55 7.1630781 -4.2148742
> 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/7xmv31353002011.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/8bzil1353002011.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/9yrp91353002011.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/10q7iu1353002011.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/11ehn71353002011.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/12ur9p1353002011.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/13rquz1353002011.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/14w0l01353002011.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/150lfk1353002011.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/16lteg1353002011.tab")
+ }
>
> try(system("convert tmp/1zumd1353002011.ps tmp/1zumd1353002011.png",intern=TRUE))
character(0)
> try(system("convert tmp/2s1cq1353002011.ps tmp/2s1cq1353002011.png",intern=TRUE))
character(0)
> try(system("convert tmp/3fx0j1353002011.ps tmp/3fx0j1353002011.png",intern=TRUE))
character(0)
> try(system("convert tmp/4378u1353002011.ps tmp/4378u1353002011.png",intern=TRUE))
character(0)
> try(system("convert tmp/5sls31353002011.ps tmp/5sls31353002011.png",intern=TRUE))
character(0)
> try(system("convert tmp/6065q1353002011.ps tmp/6065q1353002011.png",intern=TRUE))
character(0)
> try(system("convert tmp/7xmv31353002011.ps tmp/7xmv31353002011.png",intern=TRUE))
character(0)
> try(system("convert tmp/8bzil1353002011.ps tmp/8bzil1353002011.png",intern=TRUE))
character(0)
> try(system("convert tmp/9yrp91353002011.ps tmp/9yrp91353002011.png",intern=TRUE))
character(0)
> try(system("convert tmp/10q7iu1353002011.ps tmp/10q7iu1353002011.png",intern=TRUE))
character(0)
>
>
> proc.time()
user system elapsed
6.658 1.196 7.848