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
use crate::linalg::{LinalgCholesky as _, LinalgInverse as _, UPLO};
use crate::var::VarResult;
use crate::{GreenersError, VAR};
use ndarray::{s, Array2, Array3};
use std::fmt;
/// Identification scheme for SVAR.
#[derive(Debug, Clone)]
pub enum SVarIdentification {
/// Cholesky decomposition (recursive ordering).
Cholesky,
/// Short-run restrictions: A * u_t = B * e_t.
/// Masks are k x k matrices where NaN = free parameter, finite = restricted value.
ShortRun(Array2<f64>, Array2<f64>),
/// Long-run restrictions (Blanchard-Quah style).
/// Mask is k x k where NaN = free, 0.0 = restricted to zero.
LongRun(Array2<f64>),
}
/// SVAR estimation result.
#[derive(Debug)]
pub struct SVarResult {
pub var_result: VarResult,
pub a_matrix: Array2<f64>,
pub b_matrix: Array2<f64>,
pub identification: String,
}
impl SVarResult {
/// Structural impulse response function.
///
/// Returns Array3 of dimension (steps x k x k).
/// Element [h, i, j] = response of variable i to structural shock j at horizon h.
pub fn structural_irf(&self, steps: usize) -> Result<Array3<f64>, GreenersError> {
let k = self.var_result.n_vars;
let p = self.var_result.lags;
// Structural impact matrix: A^{-1} * B
let a_inv = self
.a_matrix
.inv()
.map_err(|_| GreenersError::SingularMatrix)?;
let impact = a_inv.dot(&self.b_matrix);
// Extract companion form A matrices
let mut a_matrices = Vec::new();
for l in 0..p {
let start_row = 1 + l * k;
let end_row = 1 + (l + 1) * k;
let a_lag = self
.var_result
.params
.slice(s![start_row..end_row, ..])
.t()
.to_owned();
a_matrices.push(a_lag);
}
// Compute IRF recursively
let mut phi_history = Vec::with_capacity(steps);
let mut irf = Array3::<f64>::zeros((steps, k, k));
let phi_0 = Array2::<f64>::eye(k);
let theta_0 = phi_0.dot(&impact);
irf.slice_mut(s![0, .., ..]).assign(&theta_0);
phi_history.push(phi_0);
for h in 1..steps {
let mut phi_h = Array2::<f64>::zeros((k, k));
for j in 1..=p {
if h >= j {
phi_h = phi_h + a_matrices[j - 1].dot(&phi_history[h - j]);
}
}
let theta_h = phi_h.dot(&impact);
irf.slice_mut(s![h, .., ..]).assign(&theta_h);
phi_history.push(phi_h);
}
Ok(irf)
}
/// Structural forecast error variance decomposition.
pub fn structural_fevd(&self, steps: usize) -> Result<Array3<f64>, GreenersError> {
let k = self.var_result.n_vars;
let irf = self.structural_irf(steps)?;
let mut fevd = Array3::<f64>::zeros((steps, k, k));
for i in 0..k {
let mut cum_mse = vec![0.0f64; k];
for h in 0..steps {
for j in 0..k {
cum_mse[j] += irf[[h, i, j]].powi(2);
}
let total: f64 = cum_mse.iter().sum();
if total > 1e-15 {
for j in 0..k {
fevd[[h, i, j]] = cum_mse[j] / total;
}
}
}
}
Ok(fevd)
}
}
impl fmt::Display for SVarResult {
fn fmt(&self, f: &mut fmt::Formatter<'_>) -> fmt::Result {
writeln!(
f,
"\n{:=^78}",
format!(" SVAR({}) — {} ", self.var_result.lags, self.identification)
)?;
writeln!(f, "{:<20} {:>10}", "No. Variables:", self.var_result.n_vars)?;
writeln!(f, "{:<20} {:>10}", "Observations:", self.var_result.n_obs)?;
writeln!(f, "\n{:-^78}", " A Matrix ")?;
for row in self.a_matrix.rows() {
write!(f, "[ ")?;
for val in row {
write!(f, "{:>10.4} ", val)?;
}
writeln!(f, "]")?;
}
writeln!(f, "\n{:-^78}", " B Matrix ")?;
for row in self.b_matrix.rows() {
write!(f, "[ ")?;
for val in row {
write!(f, "{:>10.4} ", val)?;
}
writeln!(f, "]")?;
}
writeln!(f, "{:=^78}", "")
}
}
/// Structural VAR model.
pub struct SVAR;
impl SVAR {
/// Fit a Structural VAR model.
///
/// * `data` — T x k data matrix
/// * `lags` — number of VAR lags
/// * `identification` — identification scheme
pub fn fit(
data: &Array2<f64>,
lags: usize,
identification: SVarIdentification,
) -> Result<SVarResult, GreenersError> {
let k = data.ncols();
// First, estimate reduced-form VAR
let var_result = VAR::fit(data, lags, None)?;
let (a_matrix, b_matrix, id_name) = match identification {
SVarIdentification::Cholesky => {
let p_chol = var_result
.sigma_u
.cholesky(UPLO::Lower)
.map_err(|_| GreenersError::SingularMatrix)?;
let a = Array2::<f64>::eye(k);
(a, p_chol, "Cholesky".to_string())
}
SVarIdentification::ShortRun(a_mask, b_mask) => {
// Solve A * Sigma_u * A' = B * B'
// With restrictions from masks (NaN = free, finite = fixed)
// Use iterative scoring or direct solve for just-identified case
let (a, b) =
solve_short_run_restrictions(&var_result.sigma_u, &a_mask, &b_mask, k)?;
(a, b, "Short-run restrictions".to_string())
}
SVarIdentification::LongRun(_c_mask) => {
// Blanchard-Quah: long-run impact matrix C(1)*A^{-1} has restrictions
// C(1) = (I - A1 - A2 - ... - Ap)^{-1}
let mut c1_inv = Array2::<f64>::eye(k);
let p = var_result.lags;
for l in 0..p {
let start_row = 1 + l * k;
let end_row = 1 + (l + 1) * k;
let a_lag = var_result
.params
.slice(s![start_row..end_row, ..])
.t()
.to_owned();
c1_inv -= &a_lag;
}
let c1 = c1_inv.inv().map_err(|_| GreenersError::SingularMatrix)?;
// Long-run impact = C(1) * A^{-1} * B
// Under Cholesky on the long-run covariance:
// C(1) * Sigma_u * C(1)' = (C(1)*P) * (C(1)*P)'
let long_run_cov = c1.dot(&var_result.sigma_u).dot(&c1.t());
let long_run_chol = long_run_cov
.cholesky(UPLO::Lower)
.map_err(|_| GreenersError::SingularMatrix)?;
// Impact matrix: P = C(1)^{-1} * long_run_chol
let b = c1_inv.dot(&long_run_chol);
let a = Array2::<f64>::eye(k);
(a, b, "Long-run restrictions".to_string())
}
};
Ok(SVarResult {
var_result,
a_matrix,
b_matrix,
identification: id_name,
})
}
}
/// Solve short-run SVAR restrictions via iterative method.
/// For just-identified case with A diagonal and B lower triangular,
/// this reduces to Cholesky-like decomposition.
fn solve_short_run_restrictions(
sigma_u: &Array2<f64>,
a_mask: &Array2<f64>,
b_mask: &Array2<f64>,
k: usize,
) -> Result<(Array2<f64>, Array2<f64>), GreenersError> {
// Initialize A and B from masks
let mut a = Array2::<f64>::eye(k);
let mut b = Array2::<f64>::eye(k);
// Apply fixed restrictions
for i in 0..k {
for j in 0..k {
if a_mask[[i, j]].is_finite() {
a[[i, j]] = a_mask[[i, j]];
}
if b_mask[[i, j]].is_finite() {
b[[i, j]] = b_mask[[i, j]];
}
}
}
// Iterative scoring: A * Sigma * A' = B * B'
// Simple iteration: given A, solve for B via Cholesky of A*Sigma*A'
for _iter in 0..100 {
let target = a.dot(sigma_u).dot(&a.t());
// B = cholesky(target)
let b_new = target
.cholesky(UPLO::Lower)
.map_err(|_| GreenersError::SingularMatrix)?;
// Apply B mask restrictions
for i in 0..k {
for j in 0..k {
if b_mask[[i, j]].is_finite() {
// Keep restricted value
} else {
b[[i, j]] = b_new[[i, j]];
}
}
}
// Check convergence
let diff = (&a.dot(sigma_u).dot(&a.t()) - &b.dot(&b.t()))
.mapv(|v| v.abs())
.sum();
if diff < 1e-10 {
break;
}
}
Ok((a, b))
}