Structural VAR

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

export type Matrix = number[][];
export type Result = Record<string, unknown>;

const round = (value:number) => Math.round(value * 1e12) / 1e12;
const rounded = (value:unknown):unknown => {
  if (typeof value === "number") return round(value);
  if (Array.isArray(value)) return value.map(rounded);
  if (value && typeof value === "object") {
    return Object.fromEntries(Object.entries(value).map(([key,item]) => [key,rounded(item)]));
  }
  return value;
};

function finiteMatrix(values:Matrix, minimumRows=1):{rows:Matrix;width:number} {
  if (!Array.isArray(values) || values.length < minimumRows) {
    throw new Error(`values must contain at least ${minimumRows} rows`);
  }
  if (!Array.isArray(values[0]) || values[0].length < 2) {
    throw new Error("values must contain at least two variables");
  }
  const width=values[0].length;
  const rows=values.map((row,rowIndex)=>{
    if(!Array.isArray(row)||row.length!==width)throw new Error("all rows must have the same width");
    return row.map((value,columnIndex)=>{
      if(typeof value!=="number"||!Number.isFinite(value))throw new Error(`values[${rowIndex}][${columnIndex}] must be finite`);
      return value;
    });
  });
  return {rows,width};
}

const identity=(size:number):Matrix=>Array.from({length:size},(_,row)=>Array.from({length:size},(_,column)=>row===column?1:0));
const transpose=(matrix:Matrix):Matrix=>matrix[0].map((_,column)=>matrix.map(row=>row[column]));
function matmul(left:Matrix,right:Matrix):Matrix{
  if(!left.length||!right.length||left[0].length!==right.length)throw new Error("matrix dimensions do not align");
  const rightT=transpose(right);
  return left.map(row=>rightT.map(column=>row.reduce((sum,value,index)=>sum+value*column[index],0)));
}
function matvec(matrix:Matrix,vector:number[]):number[]{
  if(!matrix.length||matrix[0].length!==vector.length)throw new Error("matrix and vector dimensions do not align");
  return matrix.map(row=>row.reduce((sum,value,index)=>sum+value*vector[index],0));
}
function maxAbsDifference(left:Matrix,right:Matrix):number{
  if(left.length!==right.length||left.some((row,index)=>row.length!==right[index].length))throw new Error("matrix dimensions do not align");
  return Math.max(...left.flatMap((row,rowIndex)=>row.map((value,column)=>Math.abs(value-right[rowIndex][column]))));
}
function qrDecomposition(matrix:Matrix,tolerance=1e-14):{q:Matrix;r:Matrix}{
  const size=matrix.length,columns=transpose(matrix),qColumns:Matrix=[],r:Matrix=Array.from({length:size},()=>Array(size).fill(0));
  for(let columnIndex=0;columnIndex<size;columnIndex++){
    let vector=[...columns[columnIndex]];
    for(let previous=0;previous<qColumns.length;previous++){
      const basis=qColumns[previous],projection=basis.reduce((sum,value,index)=>sum+value*vector[index],0);
      r[previous][columnIndex]=projection;
      vector=vector.map((value,index)=>value-projection*basis[index]);
    }
    let norm=Math.sqrt(vector.reduce((sum,value)=>sum+value*value,0));
    if(norm<=tolerance){
      vector=Array.from({length:size},(_,index)=>index===columnIndex?1:0);
      for(const basis of qColumns){
        const projection=basis.reduce((sum,value,index)=>sum+value*vector[index],0);
        vector=vector.map((value,index)=>value-projection*basis[index]);
      }
      norm=Math.sqrt(vector.reduce((sum,value)=>sum+value*value,0));
      if(norm<=tolerance)throw new Error("QR decomposition failed");
    }
    r[columnIndex][columnIndex]=norm;
    qColumns.push(vector.map(value=>value/norm));
  }
  return {q:transpose(qColumns),r};
}
export function companionMatrix(coefficients:Matrix[]):Matrix{
  if(!coefficients.length)throw new Error("at least one coefficient matrix is required");
  const {rows:first,width}=finiteMatrix(coefficients[0],2);
  if(first.length!==width)throw new Error("coefficient matrices must be square");
  const clean=[first,...coefficients.slice(1).map(matrix=>{const checked=finiteMatrix(matrix,width);if(checked.width!==width||checked.rows.length!==width)throw new Error("coefficient matrices must share a square shape");return checked.rows})];
  const size=width*clean.length,companion:Matrix=Array.from({length:size},()=>Array(size).fill(0));
  clean.forEach((matrix,lag)=>matrix.forEach((row,rowIndex)=>row.forEach((value,column)=>{companion[rowIndex][lag*width+column]=value})));
  for(let block=1;block<clean.length;block++)for(let index=0;index<width;index++)companion[block*width+index][(block-1)*width+index]=1;
  return companion;
}
export function companionSpectralRadius(coefficients:Matrix[]):number{
  let working=companionMatrix(coefficients);
  for(let iteration=0;iteration<384;iteration++){const {q,r}=qrDecomposition(working);working=matmul(r,q)}
  let radius=0,index=working.length-1;
  while(index>=0){
    if(index>0&&Math.abs(working[index][index-1])>1e-8){
      const a=working[index-1][index-1],b=working[index-1][index],c=working[index][index-1],d=working[index][index];
      const trace=a+d,determinant=a*d-b*c,discriminant=trace*trace-4*determinant;
      if(discriminant>=0){const root=Math.sqrt(discriminant);radius=Math.max(radius,Math.abs((trace+root)/2),Math.abs((trace-root)/2))}
      else radius=Math.max(radius,Math.sqrt(Math.abs(determinant)));
      index-=2;
    }else{radius=Math.max(radius,Math.abs(working[index][index]));index--}
  }
  return round(radius);
}
function stabilityDiagnostic(coefficients:Matrix[],boundary=0.90):Result{
  const radius=companionSpectralRadius(coefficients);
  return {companion_spectral_radius:radius,stability_state:radius>=1?"unstable":radius>=boundary?"near-boundary":"stable",stability_boundary:1,near_boundary_threshold:boundary};
}
function inverse(matrix:Matrix,tolerance=1e-12):Matrix{
  const size=matrix.length;
  if(!size||matrix.some(row=>row.length!==size))throw new Error("matrix must be square");
  const augmented=matrix.map((row,index)=>[...row,...identity(size)[index]]);
  for(let column=0;column<size;column++){
    let pivot=column;
    for(let row=column+1;row<size;row++)if(Math.abs(augmented[row][column])>Math.abs(augmented[pivot][column]))pivot=row;
    if(Math.abs(augmented[pivot][column])<=tolerance)throw new Error("matrix is singular");
    [augmented[column],augmented[pivot]]=[augmented[pivot],augmented[column]];
    const scale=augmented[column][column];
    augmented[column]=augmented[column].map(value=>value/scale);
    for(let row=0;row<size;row++){
      if(row===column)continue;
      const factor=augmented[row][column];
      augmented[row]=augmented[row].map((value,index)=>value-factor*augmented[column][index]);
    }
  }
  return augmented.map(row=>row.slice(size));
}
const ols=(design:Matrix,targets:Matrix):Matrix=>{
  const designT=transpose(design);
  return matmul(inverse(matmul(designT,design)),matmul(designT,targets));
};
function covarianceMLE(residuals:Matrix):Matrix{
  if(!residuals.length)throw new Error("residuals cannot be empty");
  const width=residuals[0].length,covariance=Array.from({length:width},()=>Array(width).fill(0));
  for(const row of residuals)for(let left=0;left<width;left++)for(let right=0;right<width;right++)covariance[left][right]+=row[left]*row[right];
  return covariance.map(row=>row.map(value=>value/residuals.length));
}

export function fitVAR(values:Matrix,lags=1,includeIntercept=true):Result{
  const {rows,width}=finiteMatrix(values,8);
  if(!Number.isInteger(lags)||lags<1||lags>=rows.length-2)throw new Error("lags must be an integer between 1 and n-3");
  const design:Matrix=[],targets:Matrix=[];
  for(let time=lags;time<rows.length;time++){
    const features:number[]=includeIntercept?[1]:[];
    for(let lag=1;lag<=lags;lag++)features.push(...rows[time-lag]);
    design.push(features);targets.push([...rows[time]]);
  }
  const beta=ols(design,targets),intercept=includeIntercept?[...beta[0]]:Array(width).fill(0),offset=includeIntercept?1:0;
  const coefficients:Matrix[]=[];
  for(let lag=0;lag<lags;lag++){
    const matrix:Matrix=Array.from({length:width},()=>Array(width).fill(0));
    for(let predictor=0;predictor<width;predictor++)for(let response=0;response<width;response++)matrix[response][predictor]=beta[offset+lag*width+predictor][response];
    coefficients.push(matrix);
  }
  const fitted=matmul(design,beta);
  const residuals=targets.map((target,index)=>target.map((value,column)=>value-fitted[index][column]));
  let oneStep=[...intercept];
  for(let lag=1;lag<=lags;lag++){
    const contribution=matvec(coefficients[lag-1],rows[rows.length-lag]);
    oneStep=oneStep.map((value,index)=>value+contribution[index]);
  }
  const rowSumBound=Math.max(...Array.from({length:width},(_,row)=>coefficients.reduce((sum,matrix)=>sum+matrix[row].reduce((inner,value)=>inner+Math.abs(value),0),0)));
  return rounded({intercept,coefficients,sigma_u_mle:covarianceMLE(residuals),residuals,fitted,one_step_forecast:oneStep,effective_observations:targets.length,variables:width,lags,row_sum_stability_bound:rowSumBound,...stabilityDiagnostic(coefficients),state:"estimated",reason:"fixed-lag-ols"}) as Result;
}

export function choleskyLower(covariance:Matrix,tolerance=1e-12):Matrix{
  const {rows:matrix,width}=finiteMatrix(covariance,2);
  if(matrix.length!==width)throw new Error("covariance must be square");
  for(let row=0;row<width;row++)for(let column=0;column<width;column++)if(Math.abs(matrix[row][column]-matrix[column][row])>1e-9)throw new Error("covariance must be symmetric");
  const lower:Matrix=Array.from({length:width},()=>Array(width).fill(0));
  for(let row=0;row<width;row++)for(let column=0;column<=row;column++){
    let subtotal=0;for(let index=0;index<column;index++)subtotal+=lower[row][index]*lower[column][index];
    if(row===column){const diagonal=matrix[row][row]-subtotal;if(diagonal<=tolerance)throw new Error("covariance must be positive definite");lower[row][column]=Math.sqrt(diagonal)}
    else lower[row][column]=(matrix[row][column]-subtotal)/lower[column][column];
  }
  return rounded(lower) as Matrix;
}

export function fitRecursiveSVAR(values:Matrix,lags=1):Result{
  const reduced=fitVAR(values,lags),impact=choleskyLower(reduced.sigma_u_mle as Matrix),inverseImpact=inverse(impact);
  const structuralShocks=(reduced.residuals as Matrix).map(row=>matvec(inverseImpact,row)),reconstructed=matmul(impact,transpose(impact)),shockCovariance=covarianceMLE(structuralShocks),width=reduced.variables as number;
  return rounded({coefficients:reduced.coefficients,intercept:reduced.intercept,sigma_u_mle:reduced.sigma_u_mle,impact_matrix:impact,structural_shocks:structuralShocks,reconstructed_sigma_u:reconstructed,covariance_reconstruction_max_error:maxAbsDifference(reduced.sigma_u_mle as Matrix,reconstructed),structural_shock_covariance_mle:shockCovariance,structural_shock_identity_max_error:maxAbsDifference(shockCovariance,identity(width)),ordering:Array.from({length:width},(_,index)=>index),companion_spectral_radius:reduced.companion_spectral_radius,stability_state:reduced.stability_state,state:"identified-recursively",reason:"lower-cholesky-ordering"}) as Result;
}

export function fitVECMFixedBeta(values:Matrix,beta:number[],differenceLags=1,includeIntercept=true):Result{
  const {rows,width}=finiteMatrix(values,12);
  if(beta.length!==width||beta.some(value=>typeof value!=="number"||!Number.isFinite(value)))throw new Error("beta must contain one finite coefficient per variable");
  if(Math.abs(beta[0]-1)>1e-12)throw new Error("beta must be normalized with beta[0] = 1");
  if(!Number.isInteger(differenceLags)||differenceLags<0||differenceLags>2)throw new Error("difference_lags must be 0, 1, or 2");
  const differences=rows.slice(1).map((row,index)=>row.map((value,column)=>value-rows[index][column]));
  const start=differenceLags+1,design:Matrix=[],targets:Matrix=[],errorCorrection:number[]=[];
  for(let time=start;time<rows.length;time++){
    const ect=beta.reduce((sum,value,column)=>sum+value*rows[time-1][column],0);
    const features:number[]=includeIntercept?[1]:[];features.push(ect);
    for(let lag=1;lag<=differenceLags;lag++)features.push(...differences[time-1-lag]);
    design.push(features);targets.push([...differences[time-1]]);errorCorrection.push(ect);
  }
  const coefficients=ols(design,targets),offset=includeIntercept?1:0,intercept=includeIntercept?[...coefficients[0]]:Array(width).fill(0),alpha=[...coefficients[offset]],gamma:Matrix[]=[];
  for(let lag=0;lag<differenceLags;lag++){
    const matrix:Matrix=Array.from({length:width},()=>Array(width).fill(0));
    for(let predictor=0;predictor<width;predictor++)for(let response=0;response<width;response++)matrix[response][predictor]=coefficients[offset+1+lag*width+predictor][response];
    gamma.push(matrix);
  }
  const fitted=matmul(design,coefficients),residuals=targets.map((target,index)=>target.map((value,column)=>value-fitted[index][column]));
  const adjustmentRoot=1+beta.reduce((sum,value,index)=>sum+value*alpha[index],0);
  const halfLife=differenceLags===0&&Math.abs(adjustmentRoot)>0&&Math.abs(adjustmentRoot)<1?-Math.log(2)/Math.log(Math.abs(adjustmentRoot)):null;
  const pi=alpha.map(value=>beta.map(betaValue=>value*betaValue));
  return rounded({beta:[...beta],alpha,gamma,intercept,sigma_u_mle:covarianceMLE(residuals),error_correction:errorCorrection,pi,adjustment_root:adjustmentRoot,error_correction_loading_root:adjustmentRoot,half_life:halfLife,exact_half_life:halfLife,half_life_scope:differenceLags===0?"exact-scalar-error-correction-channel":"not-reported-short-run-gamma-present",effective_observations:targets.length,rank:1,difference_lags:differenceLags,state:"estimated-conditional-beta",reason:"rank-one-fixed-normalized-cointegration-vector"}) as Result;
}

export function movingAverageMatrices(coefficients:Matrix[],horizon:number):Matrix[]{
  if(!coefficients.length)throw new Error("at least one coefficient matrix is required");
  const {rows:first,width}=finiteMatrix(coefficients[0],2);
  if(first.length!==width)throw new Error("coefficient matrices must be square");
  const clean=[first,...coefficients.slice(1).map(matrix=>{const checked=finiteMatrix(matrix,width);if(checked.width!==width||checked.rows.length!==width)throw new Error("coefficient matrices must share a square shape");return checked.rows})];
  if(!Number.isInteger(horizon)||horizon<0||horizon>100)throw new Error("horizon must be an integer from 0 through 100");
  const responses:Matrix[]=[identity(width)];
  for(let step=1;step<=horizon;step++){
    let current:Matrix=Array.from({length:width},()=>Array(width).fill(0));
    for(let lag=1;lag<=Math.min(clean.length,step);lag++){
      const product=matmul(clean[lag-1],responses[step-lag]);
      current=current.map((row,rowIndex)=>row.map((value,column)=>value+product[rowIndex][column]));
    }
    responses.push(current);
  }
  return rounded(responses) as Matrix[];
}

export function impulseResponses(coefficients:Matrix[],horizon:number,impactMatrix?:Matrix):Result{
  const phi=movingAverageMatrices(coefficients,horizon),width=phi[0].length,impact=impactMatrix??identity(width);
  const checked=finiteMatrix(impact,width);
  if(checked.width!==width||checked.rows.length!==width)throw new Error("impact_matrix must match coefficient width");
  const responses=phi.map(matrix=>matmul(matrix,checked.rows)),cumulative:Matrix[]=[];
  let running:Matrix=Array.from({length:width},()=>Array(width).fill(0));
  for(const matrix of responses){running=running.map((row,rowIndex)=>row.map((value,column)=>value+matrix[rowIndex][column]));cumulative.push(running.map(row=>[...row]))}
  return rounded({responses,cumulative_responses:cumulative,horizon,impact_matrix:checked.rows,...stabilityDiagnostic(coefficients),state:"computed",reason:"ma-recursion-times-declared-impact"}) as Result;
}

export function forecastErrorVarianceDecomposition(coefficients:Matrix[],sigmaU:Matrix,horizon:number,impactMatrix?:Matrix):Result{
  if(!Number.isInteger(horizon)||horizon<1||horizon>100)throw new Error("horizon must be an integer from 1 through 100");
  const impact=impactMatrix??choleskyLower(sigmaU),irf=impulseResponses(coefficients,horizon-1,impact),responses=irf.responses as Matrix[],width=responses[0].length;
  const cumulative:Matrix=Array.from({length:width},()=>Array(width).fill(0)),horizons:Matrix[]=[];
  for(const matrix of responses){
    for(let response=0;response<width;response++)for(let shock=0;shock<width;shock++)cumulative[response][shock]+=matrix[response][shock]**2;
    horizons.push(cumulative.map(row=>{const total=row.reduce((sum,value)=>sum+value,0);if(total<=0)throw new Error("forecast-error variance is zero");return row.map(value=>value/total)}));
  }
  return rounded({shares:horizons,impact_matrix:impact,horizon,row_sums:horizons.map(matrix=>matrix.map(row=>row.reduce((sum,value)=>sum+value,0))),companion_spectral_radius:irf.companion_spectral_radius,stability_state:irf.stability_state,stability_boundary:irf.stability_boundary,near_boundary_threshold:irf.near_boundary_threshold,state:"computed",reason:"orthogonalized-squared-response-share"}) as Result;
}