-
Notifications
You must be signed in to change notification settings - Fork 5
/
01classic-p4-m2.R
388 lines (341 loc) · 15.1 KB
/
01classic-p4-m2.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
rm(list = ls())
options(echo=TRUE)
# Estimation of heteroskedastic model with a lower-triangular matrix A_0 being estimated
source("./00codes/B-SVAR-MSH-identified/B-SVAR-MSH-identified-source.R")
source("./00codes/EM-MSVAR/Caspur-MSVAR-source-EM.R")
load("dataBI2015.RData")
set.seed(1234567)
Y=100*Y
p = 4
N = dim(Y)[2]
T = dim(Y)[1]
TT = T - p
M = 2
model = "classic"
# 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)
# For Transition probabilities
restrictions$M = diag(M^2)
restrictions$dj = seq(from=M, to=M, length.out=M)
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,
omega.a = 1,
omega.b = 3,
hyper.a = 1,
hyper.b = 1,
w = matrix(9 * diag(M) + matrix(1,M,M),ncol=1) # Priors: Hidden Markov Chain, transition probabilities
)
C = 0.1
# starting.values = B.SVAR.MSH.VW.identified.initialization(Y, p, M, restrictions, which.Sigma=1)
# t.0 = proc.time()
# qqq = B.SVAR.MSH.identified.Gibbs(S=100, 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.MSH.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.MSH.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.MSH.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.MSH.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.MSH.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.MSH.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=""))
# load(paste("./00estimation/classic-p",p,"-M",M,"-03.RData", sep=""))
# rm(qqq5)
# t.0 = proc.time()
# qqq1 = B.SVAR.MSH.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.MSH.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.MSH.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.MSH.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.MSH.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.MSH.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.MSH.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.MSH.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.MSH.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.MSH.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.MSH.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.MSH.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=""))
#
#
# load(paste("./00estimation/classic-p",p,"-M",M,"-09.RData", sep=""))
# rm(qqq5)
# t.0 = proc.time()
# qqq1 = B.SVAR.MSH.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.MSH.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.MSH.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.MSH.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.MSH.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.MSH.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.MSH.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.MSH.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.MSH.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.MSH.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.MSH.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.MSH.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.MSH.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.MSH.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.MSH.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.MSH.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.MSH.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.MSH.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.MSH.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.MSH.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.MSH.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.MSH.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.MSH.identified.merge(B.SVAR.MSH.identified.short.mcmc(qqq1,keep.every=10),B.SVAR.MSH.identified.short.mcmc(qqq2,keep.every=10))
rm(qqq1,qqq2)
load(paste("./00estimation/",model,"-p",p,"-M",M,"-14.RData", sep=""))
qqq = B.SVAR.MSH.identified.merge(qqq,B.SVAR.MSH.identified.short.mcmc(qqq3,keep.every=10))
qqq = B.SVAR.MSH.identified.merge(qqq,B.SVAR.MSH.identified.short.mcmc(qqq4,keep.every=10))
rm(qqq3,qqq4)
load(paste("./00estimation/",model,"-p",p,"-M",M,"-15.RData", sep=""))
qqq = B.SVAR.MSH.identified.merge(qqq,B.SVAR.MSH.identified.short.mcmc(qqq5,keep.every=10))
qqq = B.SVAR.MSH.identified.merge(qqq,B.SVAR.MSH.identified.short.mcmc(qqq6,keep.every=10))
rm(qqq5,qqq6)
t0 = proc.time()
log.mdd = B.SVAR.MSH.identified.MDD(Sp=100000,Gibbs.output=qqq, priors, restrictions)
t1 = proc.time()
t.for.time = (t1-t0)[3]
log.sddr = B.SVAR.MSH.identified.SDDR(qqq=qqq,priors)
save(log.sddr,log.mdd,qqq,priors,restrictions, file = paste("./00estimation/",model,"-p",p,"-M",M,".RData", sep=""))