外观
动态面板 GMM
差分 / 系统 GMM，语义对齐 Stata xtabond2。
核心代码
py
def _estimate(self) -> DynamicPanelGMMResult:
"""执行核心一阶/二阶 GMM 计算。"""
matrices = self._prepare_estimation_matrices()
self._prepared['estimation'] = matrices
z = matrices['z']
x = matrices['x']
y = matrices['y']
rows = matrices['rows_per_group']
instrument_count = z.shape[2]
parameter_count = x.shape[2]
h_value = _safe_int(self.dynamic_options.get('h'), 3)
h_matrix = _initial_h(
h_value,
_as_bool(self.dynamic_options.get('orthogonal'), False),
_as_bool(self.dynamic_options.get('orthogonal'), False),
matrices['system_gmm'],
self._prepared['T'],
)
zx = matrices['z_flat'].T @ matrices['x_flat']
zy = matrices['z_flat'].T @ matrices['y_flat']
first_moment_cov = np.zeros((instrument_count, instrument_count), dtype=float)
for entity in range(self._prepared['N']):
zi = z[entity]
if matrices['weight_type']:
square_root_weight = np.sqrt(matrices['weights'][entity])
weighted_zi = square_root_weight[:, None] * zi
first_moment_cov += weighted_zi.T @ h_matrix @ weighted_zi
else:
first_moment_cov += zi.T @ h_matrix @ zi
a1_unscaled, a1_rank = _symmetric_inverse(first_moment_cov)
v1_base, _inv_a1 = _symmetric_inverse(zx.T @ a1_unscaled @ zx)
b1 = v1_base @ zx.T @ a1_unscaled @ zy
residual1 = y - np.einsum('ntk,k->nt', x, b1)
sigma2 = self._compute_sigma2(matrices, residual1)
if not np.isfinite(sigma2) or sigma2 <= 0:
raise ValueError(_('估计得到的误差方差不为正，模型可能退化'))
a1 = a1_unscaled / sigma2
v1 = v1_base * sigma2
ze1 = matrices['z_flat'].T @ residual1.reshape(-1)
overid_df = max(a1_rank - parameter_count, 0)
sargan = float(ze1.T @ a1 @ ze1)
sargan_p = float(chi2.sf(sargan, overid_df)) if overid_df > 0 else np.nan
twostep = _as_bool(self.dynamic_options.get('twostep'), False)
robust = (
_as_bool(self.dynamic_options.get('robust'), False)
or bool(matrices['cluster_specs'])
or (matrices['weight_type'] == 'pweight' and not twostep)
)
one_step_nonrobust = not twostep and not robust
small = _as_bool(self.dynamic_options.get('small'), False)
if matrices['cluster_specs']:
moment_cov2, cluster_moments = self._cluster_moment_meat(
matrices['z_flat'], residual1.reshape(-1), matrices['cluster_specs'], small
)
panel_moments: List[np.ndarray] = []
else:
moment_cov2, panel_moments = self._panel_moment_meat(z, residual1)
cluster_moments = []
v1_robust = v1 @ zx.T @ a1 @ moment_cov2 @ a1 @ zx @ v1
a2, _inv_cov2 = _symmetric_inverse(moment_cov2)
v2, _inv_a2 = _symmetric_inverse(zx.T @ a2 @ zx)
b2 = v2 @ zx.T @ a2 @ zy
residual2 = y - np.einsum('ntk,k->nt', x, b2)
ze2 = matrices['z_flat'].T @ residual2.reshape(-1)
hansen = float(ze2.T @ a2 @ ze2)
hansen_p = float(chi2.sf(hansen, overid_df)) if overid_df > 0 else np.nan
v2_robust = v2.copy()
if twostep and robust:
d_matrix = self._windmeijer_d_matrix(
matrices=matrices,
a2=a2,
a2_ze=a2 @ ze2,
residual1=residual1,
panel_moments=panel_moments,
cluster_moments=cluster_moments,
small=small,
)
correction = v2 @ zx.T @ a2 @ d_matrix
v2_robust = v2 + correction @ v1_robust @ correction.T + 2.0 * correction @ v2
if twostep:
coefficients = b2
residuals = residual2
covariance = v2_robust if robust else v2
chosen_a = a2
ze = ze2
sigma2 = self._compute_sigma2(matrices, residual2)
nonrobust_covariance = v2
else:
coefficients = b1
residuals = residual1
covariance = v1_robust if robust else v1
chosen_a = a1
ze = ze1
nonrobust_covariance = v1
covariance = (covariance + covariance.T) / 2.0
ar_covariance = covariance.copy()
if small:
if one_step_nonrobust:
scale = matrices['weight_total_for_sigma'] / max(
matrices['weight_total_for_sigma'] - parameter_count, 1.0
)
elif matrices['cluster_specs']:
scale = (matrices['nobs'] - 1.0) / max(matrices['nobs'] - parameter_count, 1.0)
else:
scale = (
(matrices['nobs'] - 1.0) / max(matrices['nobs'] - parameter_count, 1.0)
* matrices['group_count'] / max(matrices['group_count'] - 1.0, 1.0)
)
covariance *= scale
sigma2 *= matrices['weight_total_for_sigma'] / max(
matrices['weight_total_for_sigma'] - parameter_count, 1.0
)
standard_errors = np.sqrt(np.clip(np.diag(covariance), 0.0, None))
test_stats = np.divide(
coefficients,
standard_errors,
out=np.full_like(coefficients, np.nan),
where=standard_errors > 0,
)
constant_count = int(matrices['include_constant'])
main_slice = slice(self._prepared['T'], 2 * self._prepared['T']) if matrices['system_gmm'] else slice(None)
main_x = x[:, main_slice, :].reshape(-1, parameter_count)
main_sample = matrices['sample'][:, main_slice].reshape(-1)
model_df = max(int(np.linalg.matrix_rank(main_x[main_sample])) - constant_count, 0)
if small:
residual_df = (
(matrices['effective_nobs'] - model_df if one_step_nonrobust else matrices['min_clusters'])
- constant_count
)
pvalues = 2.0 * student_t.sf(np.abs(test_stats), max(residual_df, 1))
else:
residual_df = np.nan
pvalues = 2.0 * norm.sf(np.abs(test_stats))
if model_df:
wald_value = float(coefficients.T @ np.linalg.pinv(covariance, hermitian=True) @ coefficients)
if small:
model_test = wald_value / model_df
model_p = float(f_dist.sf(model_test, model_df, max(residual_df, 1)))
else:
model_test = wald_value
model_p = float(chi2.sf(model_test, model_df))
else:
model_test = np.nan
model_p = np.nan
m2_vzxa = -2.0 * nonrobust_covariance @ zx.T @ chosen_a
ar_results = self._ar_tests(
matrices=matrices,
coefficients=coefficients,
residuals=residuals,
sigma2=sigma2,
z=z,
zx=zx,
chosen_a=chosen_a,
covariance=ar_covariance,
nonrobust_covariance=nonrobust_covariance,
m2_vzxa=m2_vzxa,
one_step_nonrobust=one_step_nonrobust,
)
difference_tests = self._difference_tests(
matrices=matrices,
moment_cov=first_moment_cov * sigma2 if one_step_nonrobust else moment_cov2,
zy=zy,
zx=zx,
y_flat=matrices['y_flat'],
x_flat=matrices['x_flat'],
main_stat=hansen if not one_step_nonrobust else sargan,
main_df=overid_df,
pca=bool(matrices['pca_info']),
)
fitted = np.einsum('ntk,k->nt', x, coefficients)
names = matrices['x_names']
params = pd.Series(coefficients, index=names, dtype=float)
covariance_frame = pd.DataFrame(covariance, index=names, columns=names)
std_series = pd.Series(standard_errors, index=names, dtype=float)
test_series = pd.Series(test_stats, index=names, dtype=float)
pvalue_series = pd.Series(pvalues, index=names, dtype=float)
scalars: Dict[str, Any] = {
'N': matrices['nobs'],
'sargan': sargan,
'sar_df': overid_df,
'sarganp': sargan_p,
'hansen': hansen if not one_step_nonrobust else np.nan,
'hansen_df': overid_df if not one_step_nonrobust else np.nan,
'hansenp': hansen_p if not one_step_nonrobust else np.nan,
'artests': len(ar_results),
'df_m': model_df,
'df_r': residual_df,
'chi2': model_test if not small else np.nan,
'chi2p': model_p if not small else np.nan,
'F': model_test if small else np.nan,
'F_p': model_p if small else np.nan,
'sig2': sigma2,
'sigma': math.sqrt(max(sigma2, 0.0)),
'g_min': int(np.sum(matrices['pca_info']['basis_mask'])) if matrices['pca_info'] else int(np.min(matrices['observations_per_group'][matrices['observations_per_group'] > 0])),
'g_max': int(np.sum(matrices['pca_info']['basis_mask'])) if matrices['pca_info'] else int(np.max(matrices['observations_per_group'])),
'g_avg': float(matrices['nobs'] / matrices['group_count']),
'h': h_value,
'j': a1_rank,
'j0': instrument_count,
'N_g': matrices['group_count'],
}
for order, values in ar_results.items():
scalars[f'ar{order}'] = values['stat']
scalars[f'ar{order}p'] = values['pvalue']
for index, specification in enumerate(
(item for item in matrices['cluster_specs'] if item['one_variable']), start=1
):
scalars[f'Nclust{index}'] = specification['count']
scalars.update({key: value for key, value in matrices['pca_info'].items() if key in {'components', 'kmo', 'pcaR2'}})
macros = {
'predict': 'xtab2_p',
'artype': 'levels' if _as_bool(self.dynamic_options.get('arlevels'), False) else 'first differences',
'vcetype': 'Corrected' if twostep and robust else ('Robust' if robust else ''),
'twostep': 'twostep' if twostep else '',
'small': 'small' if small else '',
'esttype': 'system' if matrices['system_gmm'] else 'difference',
'pca': 'pca' if matrices['pca_info'] else '',
'transform': 'orthogonal deviations' if _as_bool(self.dynamic_options.get('orthogonal'), False) else 'first differences',
'depvar': self.y_var,
'clustvar': ' '.join(self._prepared['cluster_vars']),
'tvar': self._prepared['time_var'],
'ivar': self._prepared['entity_var'],
'cmd': 'dynamic_panel_gmm',
'version': 'STATAU dynamic-panel-gmm 1.0',
'robust': 'robust' if robust else '',
}
if matrices['weight_type']:
macros['wtype'] = matrices['weight_type']
macros['wexp'] = f"={self._prepared['weight_var']}"
for index, metadata in enumerate(matrices['iv_meta'], start=1):
macros[f'ivinsts{index}'] = ' '.join(metadata['terms'])
iv_macro_offset = len(matrices['iv_meta'])
if matrices['include_constant']:
macros[f'ivinsts{iv_macro_offset + 1}'] = '_cons'
gmm_return_groups = self._prepared['gmm_groups']
for index, group in enumerate(gmm_return_groups, start=1):
macros[f'gmminsts{index}'] = ' '.join(term.label for term in group.terms)
for index, item in enumerate(difference_tests, start=1):
macros[f'diffgroup{index}'] = item['group']
command = self._build_command(matrices)
macros['cmdline'] = command
stored_matrices: Dict[str, np.ndarray] = {
'b': coefficients.reshape(1, -1),
'V': covariance,
'A1': a1,
'Ze': (ze1 if one_step_nonrobust else ze2).reshape(-1, 1),
}
iv_equation = [metadata['equation'] for metadata in matrices['iv_meta']]
iv_passthru = [metadata['passthru'] for metadata in matrices['iv_meta']]
iv_mz = [metadata['mz'] for metadata in matrices['iv_meta']]
if matrices['include_constant']:
iv_equation.append(0)
iv_passthru.append(0)
iv_mz.append(0)
stored_matrices['ivequation'] = np.asarray(iv_equation, dtype=float).reshape(1, -1)
stored_matrices['ivpassthru'] = np.asarray(iv_passthru, dtype=float).reshape(1, -1)
stored_matrices['ivmz'] = np.asarray(iv_mz, dtype=float).reshape(1, -1)
stored_matrices['gmmequation'] = np.asarray(
[_EQUATION_CODES[group.equation] for group in gmm_return_groups], dtype=float
).reshape(1, -1)
stored_matrices['gmmpassthru'] = np.asarray(
[int(group.passthru) for group in gmm_return_groups], dtype=float
).reshape(1, -1)
stored_matrices['gmmcollapse'] = np.asarray(
[int(group.collapse) for group in gmm_return_groups], dtype=float
).reshape(1, -1)
stored_matrices['gmmorthogonal'] = np.asarray(
[int(group.orthogonal) for group in gmm_return_groups], dtype=float
).reshape(1, -1)
stored_matrices['gmmlaglimits'] = np.asarray(
[
(group.lag_min, np.nan if group.lag_max is None else group.lag_max)
for group in gmm_return_groups
],
dtype=float,
).T if gmm_return_groups else np.empty((2, 0), dtype=float)
stored_matrices['diffsargan'] = np.asarray([
[item['restricted_stat'] for item in difference_tests],
[item['difference'] for item in difference_tests],
[item['instruments_generated'] for item in difference_tests],
[item['restricted_p'] for item in difference_tests],
[item['pvalue'] for item in difference_tests],
], dtype=float) if difference_tests else np.empty((5, 0), dtype=float)
if twostep or not one_step_nonrobust:
stored_matrices['A2'] = a2
if matrices['pca_info']:
stored_matrices['eigenvalues'] = matrices['pca_info']['eigenvalues'].reshape(1, -1)
stored_matrices['eigenvectors'] = matrices['pca_info']['eigenvectors']
if _as_bool(self.dynamic_options.get('svmat'), False):
svmat_values = {
'X': matrices['x_flat'],
'Y': matrices['y_flat'].reshape(-1, 1),
'Z': matrices['z_flat'],
'H': h_matrix,
}
if matrices['weight_type']:
svmat_values['wt'] = matrices['weights'].reshape(-1, 1)
ideqt_rows = []
for entity_index, entity_value in enumerate(self._prepared['entity_values']):
for equation_index in range(2 if matrices['system_gmm'] else 1):
for time_index in range(self._prepared['T']):
try:
entity_code = float(entity_value)
except (TypeError, ValueError):
entity_code = float(entity_index + 1)
ideqt_rows.append((entity_code, equation_index, self._prepared['minimum_time'] + time_index * self._prepared['delta']))
svmat_values['ideqt'] = np.asarray(ideqt_rows, dtype=float)
for index, specification in enumerate(
(item for item in matrices['cluster_specs'] if item['one_variable']), start=1
):
svmat_values[f'clustid{index}'] = specification['ids'].reshape(-1, 1)
stored_matrices.update(svmat_values)
return DynamicPanelGMMResult(
params=params,
covariance=covariance_frame,
std_errors=std_series,
test_stats=test_series,
pvalues=pvalue_series,
nobs=matrices['nobs'],
residuals=residuals,
fitted_values=fitted,
scalars=scalars,
macros=macros,
matrices=stored_matrices,
difference_tests=difference_tests,
instrument_names=matrices['instrument_names'],
command=command,
)