-
Notifications
You must be signed in to change notification settings - Fork 5
/
01classic-p4-m1.R
448 lines (349 loc) · 16.1 KB
/
01classic-p4-m1.R
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
411
412
413
414
415
416
417
418
419
420
421
422
423
424
425
426
427
428
429
430
431
432
433
434
435
436
437
438
439
440
441
442
443
444
445
446
447
448
rm(list = ls())
options(echo=TRUE)
# Estimation of homoskedastic model with a lower-triangular matrix A_0 being estimated
source("./00codes/B-SVAR-identified/B-SVAR-identified-source.R")
load("dataBI2015.RData")
set.seed(1234567)
Y=100*Y
p = 4
N = dim(Y)[2]
T = dim(Y)[1]
TT = T - p
model = "classic"
M = 1
# lower-triangular matrix with ones on the diagonal
r = 15
restrictions = list(
Q = matrix(0, N^2, r),
q = matrix(diag(N), ncol=1)
)
restrictions$Q[2:6,1:5] = diag(5)
restrictions$Q[9:12,6:9] = diag(4)
restrictions$Q[16:18,10:12] = diag(3)
restrictions$Q[23:24,13:14] = diag(2)
restrictions$Q[30:30,15:15] = 1
# restrictions$Q
# a = rnorm(15)
# A0 = matrix(restrictions$Q%*%a + restrictions$q, ncol=N, nrow=N)
P = matrix(0, p*N, N)
P[1:N,] = diag(N)
H.tmp = rep(0,0)
for (i in 1:p) {H.tmp = c(H.tmp, rep(1/(i^2), N))}
H = diag(H.tmp)
priors = list(
beta.P = P,
beta.H = H,
lambda.a = 1,
lambda.b = 1,
hyper.a = 1,
hyper.b = 1
)
starting.values = B.SVAR.identified.initialization(Y, p, restrictions)
C = 0.2
t.0 = proc.time()
qqq = B.SVAR.identified.Gibbs(S=1000, priors, restrictions, starting.values, debug=FALSE, print.iterations=100, C=C)
t.1 = proc.time()
(t.1-t.0)/60
t.0 = proc.time()
qqq1 = B.SVAR.identified.Gibbs(S=200000, priors, restrictions, starting.values=qqq$last.draws, debug=FALSE, print.iterations=100, C=C)
t.1 = proc.time()
t.for.time = (t.1-t.0)/60
t.0 = proc.time()
qqq2 = B.SVAR.identified.Gibbs(S=200000, priors, restrictions, starting.values=qqq1$last.draws, debug=FALSE, print.iterations=100, C=C)
t.1 = proc.time()
t.for.time = (t.1-t.0)/60
save(qqq1,qqq2,priors,restrictions,starting.values,Y,t.for.time, file = paste("./00estimation/classic-p",p,"-M",M,"-01.RData", sep=""))
rm(qqq1)
t.0 = proc.time()
qqq3 = B.SVAR.identified.Gibbs(S=200000, priors, restrictions, starting.values=qqq2$last.draws, debug=FALSE, print.iterations=100, C=C)
t.1 = proc.time()
t.for.time = (t.1-t.0)/60
rm(qqq2)
t.0 = proc.time()
qqq4 = B.SVAR.identified.Gibbs(S=200000, priors, restrictions, starting.values=qqq3$last.draws, debug=FALSE, print.iterations=100, C=C)
t.1 = proc.time()
t.for.time = (t.1-t.0)/60
save(qqq3,qqq4,priors,restrictions,starting.values,Y,t.for.time, file = paste("./00estimation/classic-p",p,"-M",M,"-02.RData", sep=""))
rm(qqq3)
t.0 = proc.time()
qqq5 = B.SVAR.identified.Gibbs(S=200000, priors, restrictions, starting.values=qqq4$last.draws, debug=FALSE, print.iterations=100, C=C)
t.1 = proc.time()
t.for.time = (t.1-t.0)/60
rm(qqq4)
t.0 = proc.time()
qqq6 = B.SVAR.identified.Gibbs(S=200000, priors, restrictions, starting.values=qqq5$last.draws, debug=FALSE, print.iterations=100, C=C)
t.1 = proc.time()
t.for.time = (t.1-t.0)/60
save(qqq5,qqq6,priors,restrictions,starting.values,Y,t.for.time, file = paste("./00estimation/classic-p",p,"-M",M,"-03.RData", sep=""))
rm(qqq5)
t.0 = proc.time()
qqq1 = B.SVAR.identified.Gibbs(S=200000, priors, restrictions, starting.values=qqq6$last.draws, debug=FALSE, print.iterations=100, C=C)
t.1 = proc.time()
t.for.time = (t.1-t.0)/60
rm(qqq6)
t.0 = proc.time()
qqq2 = B.SVAR.identified.Gibbs(S=200000, priors, restrictions, starting.values=qqq1$last.draws, debug=FALSE, print.iterations=100, C=C)
t.1 = proc.time()
t.for.time = (t.1-t.0)/60
save(qqq1,qqq2,priors,restrictions,starting.values,Y,t.for.time, file = paste("./00estimation/classic-p",p,"-M",M,"-04.RData", sep=""))
rm(qqq1)
t.0 = proc.time()
qqq3 = B.SVAR.identified.Gibbs(S=200000, priors, restrictions, starting.values=qqq2$last.draws, debug=FALSE, print.iterations=100, C=C)
t.1 = proc.time()
t.for.time = (t.1-t.0)/60
rm(qqq2)
t.0 = proc.time()
qqq4 = B.SVAR.identified.Gibbs(S=200000, priors, restrictions, starting.values=qqq3$last.draws, debug=FALSE, print.iterations=100, C=C)
t.1 = proc.time()
t.for.time = (t.1-t.0)/60
save(qqq3,qqq4,priors,restrictions,starting.values,Y,t.for.time, file = paste("./00estimation/classic-p",p,"-M",M,"-05.RData", sep=""))
rm(qqq3)
t.0 = proc.time()
qqq5 = B.SVAR.identified.Gibbs(S=200000, priors, restrictions, starting.values=qqq4$last.draws, debug=FALSE, print.iterations=100, C=C)
t.1 = proc.time()
t.for.time = (t.1-t.0)/60
rm(qqq4)
t.0 = proc.time()
qqq6 = B.SVAR.identified.Gibbs(S=200000, priors, restrictions, starting.values=qqq5$last.draws, debug=FALSE, print.iterations=100, C=C)
t.1 = proc.time()
t.for.time = (t.1-t.0)/60
save(qqq5,qqq6,priors,restrictions,starting.values,Y,t.for.time, file = paste("./00estimation/classic-p",p,"-M",M,"-06.RData", sep=""))
rm(qqq5)
t.0 = proc.time()
qqq1 = B.SVAR.identified.Gibbs(S=200000, priors, restrictions, starting.values=qqq6$last.draws, debug=FALSE, print.iterations=100, C=C)
t.1 = proc.time()
t.for.time = (t.1-t.0)/60
rm(qqq6)
t.0 = proc.time()
qqq2 = B.SVAR.identified.Gibbs(S=200000, priors, restrictions, starting.values=qqq1$last.draws, debug=FALSE, print.iterations=100, C=C)
t.1 = proc.time()
t.for.time = (t.1-t.0)/60
save(qqq1,qqq2,priors,restrictions,starting.values,Y,t.for.time, file = paste("./00estimation/classic-p",p,"-M",M,"-07.RData", sep=""))
rm(qqq1)
t.0 = proc.time()
qqq3 = B.SVAR.identified.Gibbs(S=200000, priors, restrictions, starting.values=qqq2$last.draws, debug=FALSE, print.iterations=100, C=C)
t.1 = proc.time()
t.for.time = (t.1-t.0)/60
rm(qqq2)
t.0 = proc.time()
qqq4 = B.SVAR.identified.Gibbs(S=200000, priors, restrictions, starting.values=qqq3$last.draws, debug=FALSE, print.iterations=100, C=C)
t.1 = proc.time()
t.for.time = (t.1-t.0)/60
save(qqq3,qqq4,priors,restrictions,starting.values,Y,t.for.time, file = paste("./00estimation/classic-p",p,"-M",M,"-08.RData", sep=""))
rm(qqq3)
t.0 = proc.time()
qqq5 = B.SVAR.identified.Gibbs(S=200000, priors, restrictions, starting.values=qqq4$last.draws, debug=FALSE, print.iterations=100, C=C)
t.1 = proc.time()
t.for.time = (t.1-t.0)/60
rm(qqq4)
t.0 = proc.time()
qqq6 = B.SVAR.identified.Gibbs(S=200000, priors, restrictions, starting.values=qqq5$last.draws, debug=FALSE, print.iterations=100, C=C)
t.1 = proc.time()
t.for.time = (t.1-t.0)/60
save(qqq5,qqq6,priors,restrictions,starting.values,Y,t.for.time, file = paste("./00estimation/classic-p",p,"-M",M,"-09.RData", sep=""))
rm(qqq5)
t.0 = proc.time()
qqq1 = B.SVAR.identified.Gibbs(S=200000, priors, restrictions, starting.values=qqq6$last.draws, debug=FALSE, print.iterations=100, C=C)
t.1 = proc.time()
t.for.time = (t.1-t.0)/60
rm(qqq6)
t.0 = proc.time()
qqq2 = B.SVAR.identified.Gibbs(S=200000, priors, restrictions, starting.values=qqq1$last.draws, debug=FALSE, print.iterations=100, C=C)
t.1 = proc.time()
t.for.time = (t.1-t.0)/60
save(qqq1,qqq2,priors,restrictions,starting.values,Y,t.for.time, file = paste("./00estimation/classic-p",p,"-M",M,"-10.RData", sep=""))
rm(qqq1)
t.0 = proc.time()
qqq3 = B.SVAR.identified.Gibbs(S=200000, priors, restrictions, starting.values=qqq2$last.draws, debug=FALSE, print.iterations=100, C=C)
t.1 = proc.time()
t.for.time = (t.1-t.0)/60
rm(qqq2)
t.0 = proc.time()
qqq4 = B.SVAR.identified.Gibbs(S=200000, priors, restrictions, starting.values=qqq3$last.draws, debug=FALSE, print.iterations=100, C=C)
t.1 = proc.time()
t.for.time = (t.1-t.0)/60
save(qqq3,qqq4,priors,restrictions,starting.values,Y,t.for.time, file = paste("./00estimation/classic-p",p,"-M",M,"-11.RData", sep=""))
rm(qqq3)
t.0 = proc.time()
qqq5 = B.SVAR.identified.Gibbs(S=200000, priors, restrictions, starting.values=qqq4$last.draws, debug=FALSE, print.iterations=100, C=C)
t.1 = proc.time()
t.for.time = (t.1-t.0)/60
rm(qqq4)
t.0 = proc.time()
qqq6 = B.SVAR.identified.Gibbs(S=200000, priors, restrictions, starting.values=qqq5$last.draws, debug=FALSE, print.iterations=100, C=C)
t.1 = proc.time()
t.for.time = (t.1-t.0)/60
save(qqq5,qqq6,priors,restrictions,starting.values,Y,t.for.time, file = paste("./00estimation/classic-p",p,"-M",M,"-12.RData", sep=""))
rm(qqq5)
t.0 = proc.time()
qqq1 = B.SVAR.identified.Gibbs(S=200000, priors, restrictions, starting.values=qqq6$last.draws, debug=FALSE, print.iterations=100, C=C)
t.1 = proc.time()
t.for.time = (t.1-t.0)/60
rm(qqq6)
t.0 = proc.time()
qqq2 = B.SVAR.identified.Gibbs(S=200000, priors, restrictions, starting.values=qqq1$last.draws, debug=FALSE, print.iterations=100, C=C)
t.1 = proc.time()
t.for.time = (t.1-t.0)/60
save(qqq1,qqq2,priors,restrictions,starting.values,Y,t.for.time, file = paste("./00estimation/classic-p",p,"-M",M,"-13.RData", sep=""))
rm(qqq1)
t.0 = proc.time()
qqq3 = B.SVAR.identified.Gibbs(S=200000, priors, restrictions, starting.values=qqq2$last.draws, debug=FALSE, print.iterations=100, C=C)
t.1 = proc.time()
t.for.time = (t.1-t.0)/60
rm(qqq2)
t.0 = proc.time()
qqq4 = B.SVAR.identified.Gibbs(S=200000, priors, restrictions, starting.values=qqq3$last.draws, debug=FALSE, print.iterations=100, C=C)
t.1 = proc.time()
t.for.time = (t.1-t.0)/60
save(qqq3,qqq4,priors,restrictions,starting.values,Y,t.for.time, file = paste("./00estimation/classic-p",p,"-M",M,"-14.RData", sep=""))
rm(qqq3)
t.0 = proc.time()
qqq5 = B.SVAR.identified.Gibbs(S=200000, priors, restrictions, starting.values=qqq4$last.draws, debug=FALSE, print.iterations=100, C=C)
t.1 = proc.time()
t.for.time = (t.1-t.0)/60
rm(qqq4)
t.0 = proc.time()
qqq6 = B.SVAR.identified.Gibbs(S=200000, priors, restrictions, starting.values=qqq5$last.draws, debug=FALSE, print.iterations=100, C=C)
t.1 = proc.time()
t.for.time = (t.1-t.0)/60
save(qqq5,qqq6,priors,restrictions,starting.values,Y,t.for.time, file = paste("./00estimation/classic-p",p,"-M",M,"-15.RData", sep=""))
#
#
#
# rm(qqq5)
# t.0 = proc.time()
# qqq1 = B.SVAR.identified.Gibbs(S=200000, priors, restrictions, starting.values=qqq6$last.draws, debug=FALSE, print.iterations=1000)
# t.1 = proc.time()
# t.for.time = (t.1-t.0)/60
#
# rm(qqq6)
# t.0 = proc.time()
# qqq2 = B.SVAR.identified.Gibbs(S=200000, priors, restrictions, starting.values=qqq1$last.draws, debug=FALSE, print.iterations=1000)
# t.1 = proc.time()
# t.for.time = (t.1-t.0)/60
# save(qqq1,qqq2,priors,restrictions,starting.values,Y,t.for.time, file = paste("./00estimation/",model,"-p",p,"-M",M,"-16.RData", sep=""))
#
#
# rm(qqq1)
# t.0 = proc.time()
# qqq3 = B.SVAR.identified.Gibbs(S=200000, priors, restrictions, starting.values=qqq2$last.draws, debug=FALSE, print.iterations=1000)
# t.1 = proc.time()
# t.for.time = (t.1-t.0)/60
#
# rm(qqq2)
# t.0 = proc.time()
# qqq4 = B.SVAR.identified.Gibbs(S=200000, priors, restrictions, starting.values=qqq3$last.draws, debug=FALSE, print.iterations=1000)
# t.1 = proc.time()
# t.for.time = (t.1-t.0)/60
# save(qqq3,qqq4,priors,restrictions,starting.values,Y,t.for.time, file = paste("./00estimation/",model,"-p",p,"-M",M,"-17.RData", sep=""))
#
# rm(qqq3)
# t.0 = proc.time()
# qqq5 = B.SVAR.identified.Gibbs(S=200000, priors, restrictions, starting.values=qqq4$last.draws, debug=FALSE, print.iterations=1000)
# t.1 = proc.time()
# t.for.time = (t.1-t.0)/60
#
# rm(qqq4)
# t.0 = proc.time()
# qqq6 = B.SVAR.identified.Gibbs(S=200000, priors, restrictions, starting.values=qqq5$last.draws, debug=FALSE, print.iterations=1000)
# t.1 = proc.time()
# t.for.time = (t.1-t.0)/60
# save(qqq5,qqq6,priors,restrictions,starting.values,Y,t.for.time, file = paste("./00estimation/",model,"-p",p,"-M",M,"-18.RData", sep=""))
#
# #
# rm(qqq5)
# t.0 = proc.time()
# qqq7 = B.SVAR.identified.Gibbs(S=200000, priors, restrictions, starting.values=qqq6$last.draws, debug=FALSE, print.iterations=1000)
# t.1 = proc.time()
# t.for.time = (t.1-t.0)/60
#
# rm(qqq6)
# t.0 = proc.time()
# qqq8 = B.SVAR.identified.Gibbs(S=200000, priors, restrictions, starting.values=qqq7$last.draws, debug=FALSE, print.iterations=1000)
# t.1 = proc.time()
# t.for.time = (t.1-t.0)/60
# save(qqq7,qqq8,priors,restrictions,starting.values,Y,t.for.time, file = paste("./00estimation/",model,"-p",p,"-M",M,"-19.RData", sep=""))
#
#
# rm(qqq7)
# t.0 = proc.time()
# qqq9 = B.SVAR.identified.Gibbs(S=200000, priors, restrictions, starting.values=qqq8$last.draws, debug=FALSE, print.iterations=1000)
# t.1 = proc.time()
# t.for.time = (t.1-t.0)/60
#
# rm(qqq8)
# t.0 = proc.time()
# qqq10 = B.SVAR.identified.Gibbs(S=200000, priors, restrictions, starting.values=qqq9$last.draws, debug=FALSE, print.iterations=1000)
# t.1 = proc.time()
# t.for.time = (t.1-t.0)/60
# save(qqq9,qqq10,priors,restrictions,starting.values,Y,t.for.time, file = paste("./00estimation/",model,"-p",p,"-M",M,"-20.RData", sep=""))
load(paste("./00estimation/",model,"-p",p,"-M",M,"-13.RData", sep=""))
qqq = B.SVAR.identified.merge(B.SVAR.identified.short.mcmc(qqq1,keep.every=10),B.SVAR.identified.short.mcmc(qqq2,keep.every=10))
rm(qqq1,qqq2)
load(paste("./00estimation/",model,"-p",p,"-M",M,"-14.RData", sep=""))
qqq = B.SVAR.identified.merge(qqq,B.SVAR.identified.short.mcmc(qqq3,keep.every=10))
qqq = B.SVAR.identified.merge(qqq,B.SVAR.identified.short.mcmc(qqq4,keep.every=10))
rm(qqq3,qqq4)
load(paste("./00estimation/",model,"-p",p,"-M",M,"-15.RData", sep=""))
qqq = B.SVAR.identified.merge(qqq,B.SVAR.identified.short.mcmc(qqq5,keep.every=10))
qqq = B.SVAR.identified.merge(qqq,B.SVAR.identified.short.mcmc(qqq6,keep.every=10))
rm(qqq5,qqq6)
cat("\n Computing the lnMDD\n")
t0 = proc.time()
log.mdd = B.SVAR.identified.MDD(Sp=100000,Gibbs.output=qqq, priors, restrictions)
t1 = proc.time()
t.for.time = (t1-t0)[3]
save(log.mdd,qqq,priors,restrictions,t.for.time, file = paste("./00estimation/",model,"-p",p,"-M",M,".RData", sep=""))
# p=4;M=2
# model = "classic"
# load(paste(model,"-p",p,"-M",M,"-13.RData", sep=""))
# qqq = B.SVAR.identified.merge(B.SVAR.identified.short.mcmc(qqq1,keep.every=10),B.SVAR.identified.short.mcmc(qqq2,keep.every=10))
# rm(qqq1,qqq2)
# load(paste(model,"-p",p,"-M",M,"-14.RData", sep=""))
# qqq = B.SVAR.identified.merge(qqq,B.SVAR.identified.short.mcmc(qqq3,keep.every=10))
# qqq = B.SVAR.identified.merge(qqq,B.SVAR.identified.short.mcmc(qqq4,keep.every=10))
# rm(qqq3,qqq4)
# load(paste(model,"-p",p,"-M",M,"-15.RData", sep=""))
# qqq = B.SVAR.identified.merge(qqq,B.SVAR.identified.short.mcmc(qqq5,keep.every=10))
# qqq = B.SVAR.identified.merge(qqq,B.SVAR.identified.short.mcmc(qqq6,keep.every=10))
# rm(qqq5,qqq6)
#
#
# qqq = B.SVAR.identified.short.mcmc(qqq,keep.every=1)
#
#
#
# t0 = proc.time()
# log.mdd = B.SVAR.identified.MDD(Sp=100000,Gibbs.output=qqq, priors, restrictions)
# t1 = proc.time()
# t.for.time = (t1-t0)[3]
# save(log.mdd,qqq,priors,restrictions,t.for.time, file = paste("./00estimation/mdd-",model,"-p",p,"-M",M,".RData", sep=""))
# load(paste("./00estimation/mdd-",model,"-p",p,"-M",M,".RData", sep=""))
# log.sddr = B.SVAR.identified.SDDR(qqq=qqq,priors)
# save(log.sddr,qqq,priors, file = paste("./00estimation/sddr-",model,"-p",p,"-M",M,".RData", sep=""))
# head(Y)
# qqq = qqq1
# plot.ts(qqq$posteriors$MH.logkernel)
# 1-rejectionRate(as.mcmc(qqq$posteriors$a[1,]))
#
# plot.ts(t(qqq$posteriors$a[1:10,]))
# plot.ts(t(qqq$posteriors$a[11:20,]))
# plot.ts(t(qqq$posteriors$a[21:30,]))
#
# plot.ts(t(qqq$posteriors$lambda.omega[,1,]))
# plot.ts(t(qqq$posteriors$lambda.omega[,2,]))
# plot.ts(t(qqq$posteriors$lambda.omega[,3,]))
# apply(qqq$posteriors$lambda.omega,1:2,mean)
#
# plot.ts(t(qqq$posteriors$PR_TR[1,,]))
# plot.ts(t(qqq$posteriors$PR_TR[2,,]))
# plot.ts(t(qqq$posteriors$PR_TR[3,,]))
# apply(qqq$posteriors$PR_TR,1:2,mean)
#
# plot.ts(t(qqq$posteriors$hyper))
# apply(qqq$posteriors$hyper,1,mean)
#
# plot.ts(t(qqq$posteriors$alpha[1,,]))
# apply(qqq$posteriors$alpha,1:2,mean)
#
# plot.ts(G.probabilities.regime(posteriors=qqq$posteriors, keep.iterations=1:100000))