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
+ ,195
+ ,232
+ ,50
+ ,81
+ ,125
+ ,183
+ ,93
+ ,8841
+ ,201.2
+ ,206
+ ,256
+ ,59
+ ,95
+ ,98
+ ,162
+ ,68
+ ,8733
+ ,213.9
+ ,235
+ ,242
+ ,52
+ ,93
+ ,100
+ ,179
+ ,77
+ ,8885
+ ,209.7
+ ,283
+ ,302
+ ,66
+ ,109
+ ,93
+ ,162
+ ,78
+ ,8933
+ ,202.4
+ ,352
+ ,282
+ ,62
+ ,103
+ ,123
+ ,206
+ ,95
+ ,8854
+ ,187.8
+ ,358
+ ,288
+ ,59
+ ,101
+ ,116
+ ,194
+ ,88
+ ,8748
+ ,173.7
+ ,396
+ ,321
+ ,70
+ ,121
+ ,124
+ ,198
+ ,88
+ ,8827
+ ,172.3
+ ,398
+ ,316
+ ,74
+ ,112
+ ,126
+ ,219
+ ,102
+ ,8850
+ ,148
+ ,415
+ ,396
+ ,84
+ ,151
+ ,126
+ ,212
+ ,103
+ ,8761
+ ,129.8
+ ,385
+ ,362
+ ,71
+ ,142
+ ,156
+ ,265
+ ,131
+ ,8617
+ ,129.8
+ ,453
+ ,392
+ ,81
+ ,144
+ ,141
+ ,234
+ ,127
+ ,8758
+ ,117.9
+ ,555
+ ,414
+ ,92
+ ,154
+ ,163
+ ,259
+ ,133
+ ,8806
+ ,112.1
+ ,490
+ ,417
+ ,89
+ ,164
+ ,164
+ ,287
+ ,127
+ ,8710
+ ,94
+ ,554
+ ,476
+ ,100
+ ,188
+ ,156
+ ,278
+ ,138
+ ,8681
+ ,102.4
+ ,607
+ ,488
+ ,103
+ ,189
+ ,180
+ ,317
+ ,158
+ ,8819
+ ,115.8
+ ,711
+ ,489
+ ,97
+ ,188
+ ,187
+ ,320
+ ,167
+ ,8834
+ ,122.9
+ ,619
+ ,467
+ ,107
+ ,185
+ ,194
+ ,326
+ ,162
+ ,8742
+ ,120.9
+ ,744
+ ,460
+ ,93
+ ,188
+ ,168
+ ,316
+ ,149
+ ,8766
+ ,128.4
+ ,650
+ ,482
+ ,97
+ ,200
+ ,170
+ ,306
+ ,153
+ ,8902
+ ,148.8
+ ,688
+ ,510
+ ,100
+ ,211
+ ,177
+ ,315
+ ,166
+ ,8980
+ ,141.3
+ ,834
+ ,493
+ ,89
+ ,202
+ ,189
+ ,329
+ ,179
+ ,9031
+ ,163.7
+ ,882
+ ,476
+ ,102
+ ,198
+ ,194
+ ,316
+ ,176
+ ,8984
+ ,165.3
+ ,756
+ ,448
+ ,96
+ ,189
+ ,170
+ ,316
+ ,159
+ ,9150
+ ,187.3
+ ,830
+ ,410
+ ,81
+ ,174
+ ,156
+ ,297
+ ,151
+ ,9231
+ ,209.7
+ ,871
+ ,466
+ ,91
+ ,199
+ ,148
+ ,266
+ ,143
+ ,9230
+ ,230.1
+ ,821
+ ,417
+ ,84
+ ,175
+ ,167
+ ,296
+ ,169
+ ,9194
+ ,230.3
+ ,776
+ ,387
+ ,78
+ ,160
+ ,150
+ ,275
+ ,141
+ ,9307
+ ,234.9
+ ,770
+ ,370
+ ,70
+ ,160
+ ,141
+ ,252
+ ,134
+ ,9336
+ ,238.3
+ ,988
+ ,344
+ ,67
+ ,145
+ ,134
+ ,239
+ ,130
+ ,9310
+ ,246.8
+ ,873
+ ,396
+ ,76
+ ,172
+ ,127
+ ,231
+ ,112
+ ,9236
+ ,249.3
+ ,775
+ ,349
+ ,65
+ ,147
+ ,142
+ ,256
+ ,141
+ ,9244
+ ,247
+ ,712
+ ,326
+ ,66
+ ,138
+ ,132
+ ,232
+ ,116
+ ,9222
+ ,244.9
+ ,637
+ ,303
+ ,62
+ ,122
+ ,118
+ ,230
+ ,95
+ ,9186
+ ,246.7
+ ,492
+ ,300
+ ,66
+ ,118
+ ,115
+ ,205
+ ,98
+ ,9095
+ ,197.4
+ ,543
+ ,329
+ ,68
+ ,133
+ ,113
+ ,195
+ ,104
+ ,9216
+ ,153.9
+ ,540
+ ,304
+ ,59
+ ,118
+ ,123
+ ,207
+ ,121
+ ,9237
+ ,128.4
+ ,676
+ ,286
+ ,68
+ ,112
+ ,123
+ ,197
+ ,106
+ ,9207
+ ,130.7
+ ,645
+ ,281
+ ,68
+ ,109
+ ,103
+ ,194
+ ,90
+ ,9189
+ ,125.4
+ ,583
+ ,377
+ ,84
+ ,152
+ ,101
+ ,181
+ ,99
+ ,9183
+ ,115.6
+ ,615
+ ,344
+ ,75
+ ,141
+ ,135
+ ,246
+ ,130
+ ,9277
+ ,117.5
+ ,635
+ ,369
+ ,79
+ ,144
+ ,122
+ ,220
+ ,123
+ ,9305
+ ,125.3
+ ,579
+ ,390
+ ,92
+ ,152
+ ,142
+ ,234
+ ,133
+ ,9268
+ ,128.3
+ ,593
+ ,406
+ ,88
+ ,172
+ ,140
+ ,264
+ ,126
+ ,9259
+ ,134.7
+ ,547
+ ,426
+ ,98
+ ,168
+ ,138
+ ,266
+ ,137
+ ,9197
+ ,134.7
+ ,594
+ ,467
+ ,104
+ ,185
+ ,153
+ ,282
+ ,142
+ ,9293
+ ,134.1
+ ,601
+ ,437
+ ,95
+ ,174
+ ,172
+ ,312
+ ,153
+ ,9270
+ ,122.7
+ ,565
+ ,410
+ ,99
+ ,159
+ ,160
+ ,297
+ ,138
+ ,9395
+ ,117.8
+ ,631
+ ,390
+ ,93
+ ,155
+ ,146
+ ,269
+ ,139
+ ,9316
+ ,109.1
+ ,223
+ ,418
+ ,102
+ ,171
+ ,136
+ ,252
+ ,137
+ ,9237
+ ,108
+ ,234
+ ,398
+ ,91
+ ,161
+ ,139
+ ,265
+ ,152
+ ,9207
+ ,134.7
+ ,223
+ ,422
+ ,105
+ ,173
+ ,139
+ ,246
+ ,151
+ ,9189
+ ,134.7
+ ,680
+ ,439
+ ,100
+ ,179
+ ,140
+ ,263
+ ,158
+ ,9183
+ ,134.1
+ ,161
+ ,419
+ ,99
+ ,171
+ ,150
+ ,274
+ ,162
+ ,9277
+ ,122.7
+ ,234
+ ,484
+ ,111
+ ,202
+ ,142
+ ,262
+ ,156
+ ,9305
+ ,117.8
+ ,191
+ ,491
+ ,110
+ ,199
+ ,130
+ ,298
+ ,186
+ ,9268
+ ,109.1
+ ,586)
+ ,dim=c(9
+ ,56)
+ ,dimnames=list(c('werkeloosheid'
+ ,'onderwijshoog'
+ ,'onderwijsmiddelbaar'
+ ,'onderwijslaag'
+ ,'autochtoon'
+ ,'allochtonen'
+ ,'banen'
+ ,'vacatures'
+ ,'faillietevenootschappen')
+ ,1:56))
> y <- array(NA,dim=c(9,56),dimnames=list(c('werkeloosheid','onderwijshoog','onderwijsmiddelbaar','onderwijslaag','autochtoon','allochtonen','banen','vacatures','faillietevenootschappen'),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 faillietevenootschappen
1 79 8909 201.0 195
2 93 8841 201.2 206
3 68 8733 213.9 235
4 77 8885 209.7 283
5 78 8933 202.4 352
6 95 8854 187.8 358
7 88 8748 173.7 396
8 88 8827 172.3 398
9 102 8850 148.0 415
10 103 8761 129.8 385
11 131 8617 129.8 453
12 127 8758 117.9 555
13 133 8806 112.1 490
14 127 8710 94.0 554
15 138 8681 102.4 607
16 158 8819 115.8 711
17 167 8834 122.9 619
18 162 8742 120.9 744
19 149 8766 128.4 650
20 153 8902 148.8 688
21 166 8980 141.3 834
22 179 9031 163.7 882
23 176 8984 165.3 756
24 159 9150 187.3 830
25 151 9231 209.7 871
26 143 9230 230.1 821
27 169 9194 230.3 776
28 141 9307 234.9 770
29 134 9336 238.3 988
30 130 9310 246.8 873
31 112 9236 249.3 775
32 141 9244 247.0 712
33 116 9222 244.9 637
34 95 9186 246.7 492
35 98 9095 197.4 543
36 104 9216 153.9 540
37 121 9237 128.4 676
38 106 9207 130.7 645
39 90 9189 125.4 583
40 99 9183 115.6 615
41 130 9277 117.5 635
42 123 9305 125.3 579
43 133 9268 128.3 593
44 126 9259 134.7 547
45 137 9197 134.7 594
46 142 9293 134.1 601
47 153 9270 122.7 565
48 138 9395 117.8 631
49 139 9316 109.1 223
50 137 9237 108.0 234
51 152 9207 134.7 223
52 151 9189 134.7 680
53 158 9183 134.1 161
54 162 9277 122.7 234
55 156 9305 117.8 191
56 186 9268 109.1 586
> k <- length(x[1,])
> df <- as.data.frame(x)
> (mylm <- lm(df))
Call:
lm(formula = df)
Coefficients:
(Intercept) onderwijshoog onderwijsmiddelbaar
345.737282 0.779542 1.741148
onderwijslaag autochtoon allochtonen
0.083906 -0.010631 0.053791
banen vacatures faillietevenootschappen
-0.032944 -0.075656 0.003553
> (mysum <- summary(mylm))
Call:
lm(formula = df)
Residuals:
Min 1Q Median 3Q Max
-13.4777 -6.7095 -0.9929 6.2833 17.1687
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 345.737282 67.341367 5.134 5.34e-06 ***
onderwijshoog 0.779542 0.238247 3.272 0.00201 **
onderwijsmiddelbaar 1.741148 0.129466 13.449 < 2e-16 ***
onderwijslaag 0.083906 0.165126 0.508 0.61374
autochtoon -0.010631 0.105454 -0.101 0.92013
allochtonen 0.053791 0.128484 0.419 0.67737
banen -0.032944 0.007442 -4.427 5.67e-05 ***
vacatures -0.075656 0.038890 -1.945 0.05772 .
faillietevenootschappen 0.003553 0.008069 0.440 0.66176
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 8.533 on 47 degrees of freedom
Multiple R-squared: 0.9887, Adjusted R-squared: 0.9868
F-statistic: 513.7 on 8 and 47 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.6374135 0.7251730 0.36258650
[2,] 0.5728996 0.8542008 0.42710040
[3,] 0.4339297 0.8678594 0.56607030
[4,] 0.3436606 0.6873212 0.65633941
[5,] 0.4269829 0.8539657 0.57301714
[6,] 0.6002627 0.7994745 0.39973725
[7,] 0.5860851 0.8278299 0.41391493
[8,] 0.5248666 0.9502668 0.47513342
[9,] 0.6195631 0.7608738 0.38043690
[10,] 0.5844125 0.8311750 0.41558751
[11,] 0.6419110 0.7161780 0.35808902
[12,] 0.7677188 0.4645624 0.23228121
[13,] 0.7871515 0.4256970 0.21284850
[14,] 0.7509843 0.4980314 0.24901568
[15,] 0.6889309 0.6221381 0.31106906
[16,] 0.6115626 0.7768748 0.38843741
[17,] 0.5456081 0.9087838 0.45439188
[18,] 0.5191125 0.9617750 0.48088750
[19,] 0.4387534 0.8775068 0.56124658
[20,] 0.3717776 0.7435552 0.62822242
[21,] 0.2952780 0.5905559 0.70472203
[22,] 0.3312569 0.6625138 0.66874311
[23,] 0.3135070 0.6270140 0.68649298
[24,] 0.2826811 0.5653622 0.71731892
[25,] 0.3669453 0.7338905 0.63305473
[26,] 0.3485885 0.6971771 0.65141146
[27,] 0.3126087 0.6252173 0.68739133
[28,] 0.2536726 0.5073453 0.74632737
[29,] 0.2268099 0.4536198 0.77319008
[30,] 0.2200685 0.4401370 0.77993148
[31,] 0.2490644 0.4981288 0.75093559
[32,] 0.8415439 0.3169121 0.15845607
[33,] 0.9126806 0.1746388 0.08731941
> postscript(file="/var/wessaorg/rcomp/tmp/1gpsj1353341076.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/2wr6e1353341076.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/3dzbp1353341076.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/4fyjn1353341076.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/533301353341076.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
15.22066705 -1.54190863 -8.24704432 -9.25991849 13.10485488 -0.02246097
7 8 9 10 11 12
7.94091672 -0.59614677 5.11724339 5.08644715 -11.55464067 11.69430346
13 14 15 16 17 18
7.17840516 -9.11735364 -1.00416181 9.43172778 17.16873293 -11.28412038
19 20 21 22 23 24
-8.12740972 -4.73982097 2.06151572 10.95381502 -10.59168328 -8.44562151
25 26 27 28 29 30
-3.01471870 4.12127303 -1.31923830 6.27657917 -3.16177443 0.15346248
31 32 33 34 35 36
-2.30344090 0.56081634 -6.23689371 3.48800105 -2.74997716 -0.98165244
37 38 39 40 41 42
3.80492738 -10.76768717 -8.81107210 -1.66039109 -9.15175761 10.40905793
43 44 45 46 47 48
4.23692615 -10.25346561 6.30354983 14.76172658 7.57093498 8.73710336
49 50 51 52 53 54
-0.60912971 -9.44235634 -3.30552880 -13.47766625 -1.70577015 -5.95939529
55 56
-2.71856437 6.77978373
> postscript(file="/var/wessaorg/rcomp/tmp/6m3sf1353341076.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 15.22066705 NA
1 -1.54190863 15.22066705
2 -8.24704432 -1.54190863
3 -9.25991849 -8.24704432
4 13.10485488 -9.25991849
5 -0.02246097 13.10485488
6 7.94091672 -0.02246097
7 -0.59614677 7.94091672
8 5.11724339 -0.59614677
9 5.08644715 5.11724339
10 -11.55464067 5.08644715
11 11.69430346 -11.55464067
12 7.17840516 11.69430346
13 -9.11735364 7.17840516
14 -1.00416181 -9.11735364
15 9.43172778 -1.00416181
16 17.16873293 9.43172778
17 -11.28412038 17.16873293
18 -8.12740972 -11.28412038
19 -4.73982097 -8.12740972
20 2.06151572 -4.73982097
21 10.95381502 2.06151572
22 -10.59168328 10.95381502
23 -8.44562151 -10.59168328
24 -3.01471870 -8.44562151
25 4.12127303 -3.01471870
26 -1.31923830 4.12127303
27 6.27657917 -1.31923830
28 -3.16177443 6.27657917
29 0.15346248 -3.16177443
30 -2.30344090 0.15346248
31 0.56081634 -2.30344090
32 -6.23689371 0.56081634
33 3.48800105 -6.23689371
34 -2.74997716 3.48800105
35 -0.98165244 -2.74997716
36 3.80492738 -0.98165244
37 -10.76768717 3.80492738
38 -8.81107210 -10.76768717
39 -1.66039109 -8.81107210
40 -9.15175761 -1.66039109
41 10.40905793 -9.15175761
42 4.23692615 10.40905793
43 -10.25346561 4.23692615
44 6.30354983 -10.25346561
45 14.76172658 6.30354983
46 7.57093498 14.76172658
47 8.73710336 7.57093498
48 -0.60912971 8.73710336
49 -9.44235634 -0.60912971
50 -3.30552880 -9.44235634
51 -13.47766625 -3.30552880
52 -1.70577015 -13.47766625
53 -5.95939529 -1.70577015
54 -2.71856437 -5.95939529
55 6.77978373 -2.71856437
56 NA 6.77978373
> dum1 <- dum[2:length(myerror),]
> dum1
lag(myerror, k = 1) myerror
[1,] -1.54190863 15.22066705
[2,] -8.24704432 -1.54190863
[3,] -9.25991849 -8.24704432
[4,] 13.10485488 -9.25991849
[5,] -0.02246097 13.10485488
[6,] 7.94091672 -0.02246097
[7,] -0.59614677 7.94091672
[8,] 5.11724339 -0.59614677
[9,] 5.08644715 5.11724339
[10,] -11.55464067 5.08644715
[11,] 11.69430346 -11.55464067
[12,] 7.17840516 11.69430346
[13,] -9.11735364 7.17840516
[14,] -1.00416181 -9.11735364
[15,] 9.43172778 -1.00416181
[16,] 17.16873293 9.43172778
[17,] -11.28412038 17.16873293
[18,] -8.12740972 -11.28412038
[19,] -4.73982097 -8.12740972
[20,] 2.06151572 -4.73982097
[21,] 10.95381502 2.06151572
[22,] -10.59168328 10.95381502
[23,] -8.44562151 -10.59168328
[24,] -3.01471870 -8.44562151
[25,] 4.12127303 -3.01471870
[26,] -1.31923830 4.12127303
[27,] 6.27657917 -1.31923830
[28,] -3.16177443 6.27657917
[29,] 0.15346248 -3.16177443
[30,] -2.30344090 0.15346248
[31,] 0.56081634 -2.30344090
[32,] -6.23689371 0.56081634
[33,] 3.48800105 -6.23689371
[34,] -2.74997716 3.48800105
[35,] -0.98165244 -2.74997716
[36,] 3.80492738 -0.98165244
[37,] -10.76768717 3.80492738
[38,] -8.81107210 -10.76768717
[39,] -1.66039109 -8.81107210
[40,] -9.15175761 -1.66039109
[41,] 10.40905793 -9.15175761
[42,] 4.23692615 10.40905793
[43,] -10.25346561 4.23692615
[44,] 6.30354983 -10.25346561
[45,] 14.76172658 6.30354983
[46,] 7.57093498 14.76172658
[47,] 8.73710336 7.57093498
[48,] -0.60912971 8.73710336
[49,] -9.44235634 -0.60912971
[50,] -3.30552880 -9.44235634
[51,] -13.47766625 -3.30552880
[52,] -1.70577015 -13.47766625
[53,] -5.95939529 -1.70577015
[54,] -2.71856437 -5.95939529
[55,] 6.77978373 -2.71856437
> z <- as.data.frame(dum1)
> z
lag(myerror, k = 1) myerror
1 -1.54190863 15.22066705
2 -8.24704432 -1.54190863
3 -9.25991849 -8.24704432
4 13.10485488 -9.25991849
5 -0.02246097 13.10485488
6 7.94091672 -0.02246097
7 -0.59614677 7.94091672
8 5.11724339 -0.59614677
9 5.08644715 5.11724339
10 -11.55464067 5.08644715
11 11.69430346 -11.55464067
12 7.17840516 11.69430346
13 -9.11735364 7.17840516
14 -1.00416181 -9.11735364
15 9.43172778 -1.00416181
16 17.16873293 9.43172778
17 -11.28412038 17.16873293
18 -8.12740972 -11.28412038
19 -4.73982097 -8.12740972
20 2.06151572 -4.73982097
21 10.95381502 2.06151572
22 -10.59168328 10.95381502
23 -8.44562151 -10.59168328
24 -3.01471870 -8.44562151
25 4.12127303 -3.01471870
26 -1.31923830 4.12127303
27 6.27657917 -1.31923830
28 -3.16177443 6.27657917
29 0.15346248 -3.16177443
30 -2.30344090 0.15346248
31 0.56081634 -2.30344090
32 -6.23689371 0.56081634
33 3.48800105 -6.23689371
34 -2.74997716 3.48800105
35 -0.98165244 -2.74997716
36 3.80492738 -0.98165244
37 -10.76768717 3.80492738
38 -8.81107210 -10.76768717
39 -1.66039109 -8.81107210
40 -9.15175761 -1.66039109
41 10.40905793 -9.15175761
42 4.23692615 10.40905793
43 -10.25346561 4.23692615
44 6.30354983 -10.25346561
45 14.76172658 6.30354983
46 7.57093498 14.76172658
47 8.73710336 7.57093498
48 -0.60912971 8.73710336
49 -9.44235634 -0.60912971
50 -3.30552880 -9.44235634
51 -13.47766625 -3.30552880
52 -1.70577015 -13.47766625
53 -5.95939529 -1.70577015
54 -2.71856437 -5.95939529
55 6.77978373 -2.71856437
> 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/7419r1353341076.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/8bye91353341076.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/9j2xz1353341076.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/10v0px1353341076.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/11b3ji1353341076.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/1260fw1353341076.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/13f58r1353341076.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/14i1rp1353341076.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/157o5n1353341076.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/16foxr1353341076.tab")
+ }
>
> try(system("convert tmp/1gpsj1353341076.ps tmp/1gpsj1353341076.png",intern=TRUE))
character(0)
> try(system("convert tmp/2wr6e1353341076.ps tmp/2wr6e1353341076.png",intern=TRUE))
character(0)
> try(system("convert tmp/3dzbp1353341076.ps tmp/3dzbp1353341076.png",intern=TRUE))
character(0)
> try(system("convert tmp/4fyjn1353341076.ps tmp/4fyjn1353341076.png",intern=TRUE))
character(0)
> try(system("convert tmp/533301353341076.ps tmp/533301353341076.png",intern=TRUE))
character(0)
> try(system("convert tmp/6m3sf1353341076.ps tmp/6m3sf1353341076.png",intern=TRUE))
character(0)
> try(system("convert tmp/7419r1353341076.ps tmp/7419r1353341076.png",intern=TRUE))
character(0)
> try(system("convert tmp/8bye91353341076.ps tmp/8bye91353341076.png",intern=TRUE))
character(0)
> try(system("convert tmp/9j2xz1353341076.ps tmp/9j2xz1353341076.png",intern=TRUE))
character(0)
> try(system("convert tmp/10v0px1353341076.ps tmp/10v0px1353341076.png",intern=TRUE))
character(0)
>
>
> proc.time()
user system elapsed
6.329 0.878 7.225