Rolling variance estimator

Dear all,

Currently I am trying to implement a rolling estimator for volatility, and simplest way of doing this(before moving on to more complicated measures) was to implement a rolling variance estimator, based on exponential weighting.

I have coded this example from “Incremental calculation of weighted mean and variance” www-uxsup.csx.cam.ac.uk/~fanf2/hermes/doc/antiforgery/stats.pdf into dynare as follows:

ewma = (1-alph)*ewma(-1) + alph*y; ewmv = (1-alph)*ewmv(-1) + alph*((y-ewma)*(y-ewma(-1)));

When I run this in dynare, I do not get a value for the variance estimator but I do get the moving average estimate (and its impulse response function).
The variance is 0 everywhere according to Dynare, while the variance of y is reported to be 127.

When I exported the raw data using dynasave, and analysed it in mathematica, using the exact same estimator, it did work.

So would it be possible to calculate variance of a variable while the simulation is running?

With regards to the rolling variance estimator,

I had tried to implement it in vermandel 2013 canonical NK model.
Where the aggressiveness of monetary policy is dependant on the amount of volatility in the simulated economy.
The reasoning behind this is since volatility is a sign of destabilisation, the central bank would act differently during the shock vs in the steady state.
Currently, this is based on the volatility in the output gap (y).

So to do this, I had set up a volatility estimator(based on the stats pdf and extended this to use the correlation function between slightly different time windows). Simple one pass estimator works in two steps, 1) calculate the incrementally changed exponentially weighted moving average, 2) calculate the incrementally changed exponentially weighted moving variance.

[code]

  1. ewma = (1-alph)ewma(-1) + alphy;
  2. ewmv = (1-alph)ewmv(-1) + alph((y-ewma)*(y-ewma(-1)));[/code]

This however doesn’t work in dynare, also the estimator I want to use is slightly more involved, since it takes to moving averages(eq 3 and 4) , then calculates two incrementally updated variances( eq. 5 and 6 ) and then uses eq. 3 and eq.4 to calculate incrementally updated covariance (eq. 7) between slightly different time windows and then combines eq 7 and eq 5 and eq 6 into one estimate for the correlation (eq. 8):

Rolling volatility estimator, where the alph_'s are predetermined

3)        yema = (1-alph_a)*yema(-1) + alph_a*y;
4)        yemb = (1-alph_b)*yemb(-1) + alph_b*y;
5)        ysema = (1-alph_a)*ysema(-1) + alph_a*((y-yema)*(y-yema(-1)));
6)        ysemb = (1-alph_b)*ysemb(-1) + alph_b*((y-yemb)*(y-yemb(-1)));
7)        ycove = (1-alph_c)*(ycove(-1) + alph_c*((y-yema(-1))*(y-yemb(-1))));
8)        ycore = ycove/(sqrt(ysemb)*sqrt(ysema))

Equations 3 and 4 work, but 5 to 8 do not, where in matlab or mathematica, they do and produce a very clear graph.
volatility plot.pdf (22.2 KB)
Then I use eq 8 as a measure of volatility, and plug this into a calibrated logistic function (eq. 9) outputs gives a low (predetermined) monetary policy coefficient(which is denoted by phi_dy in eq 10) value of (cay - delt) when there is low volatility(chi_min), and a high coefficent (cay + delt) when measured volatility is high (when ycore=chi_max).
So I can vary the “central” monetary policy coefficient value(which is cay) and vary how aggressive monetary policy would react under either high/low volatility, by varying delta. So by looping over varying amounts of delta I should be able to get some sense whether this approach would help.

9)   phi_y = ((-cay - delt)*(exp(kay*ycore + (kay*(chi_max))/2)) + 
      (cay - delt)*(exp(kay*ycore + (kay*(chi_min))/2)) - 
      (cay - delt)*(exp(kay*(chi_max) + (kay*(chi_min))/2)) + 
      (cay + delt)*(exp((kay*(chi_max))/2 + kay*(chi_min))))/
      ((-(exp((kay*(chi_max))/2)) + exp((kay*(chi_min))/2))*
      (exp(kay*ycore) + exp((1/2)*kay*((chi_max) + (chi_min)))));

(This code has been derived and checked using mathematica)

And then the final piece is putting this to work in the monetary policy, which is simple OSR. So an Interest rule(eq. 10) which react to inflation and output gap, and some smoothing parameter. where right now, only the output gap responds to its volatility.

10)    r = rho*r(-1) +  (1-rho)*(phi_pi*pi + phi_y*y) + s_r  ;

The rest of the economy is characterized by a IS curve (eq. 11), which determines the output gap, and an AS curve (eq.12), which determines inflation.

11) y = y(+1) - 1/sigmaC*(r-pi(+1)) + s_b; 12) pi = beta*pi(+1) + ((1-theta)*(1-beta*theta)/theta)*(sigmaC+sigmaL)*y + s_p;
The model is simulated under 3 shocks,

s_b = rho_b*s_b(-1)+e_b; s_p = rho_p*s_p(-1)+e_p-u*e_p(-1); s_r = rho_r*s_r(-1)+e_r; shocks; var e_b; stderr 1; var e_p; stderr 1; var e_r; stderr 1; end;

I have supplied dynare with some initial starting values for the estimators (although that should be unnecessary, since it shouldn’t look for a steady state in there).

initval; yema = 0.6; yemb = 0.5; ysema = 1; ysemb = 1; ycove = 0.5; ycore = 0.5; phi_y = 0.2; end;

The model is simulated with stochastic simulation in third order and followed by the dynasave command as such:
stoch_simul(order=3,irf=30,periods=200) y r pi phi_y ycore;
dynasave(ma);

After running this, the volatility estimator(ycore) gives no results, therefore phi_y doesnt respond to it either. I had hoped since it tried to find some steady state, completing the code would work, but it doesn’t:-(…
When I load the saved data, I can see the estimator works, but it would have been really nice to actually see the decomposed impulse response to volatility.
Does anyone know why this(the rolling estimation of volatility) doesn’t work?
In the mean time, I am trying to implement this in a different paper with observation equations (Uribe 2007 extended), and I am hoping this will yield results, since if volatility is observed, there should be no reason for dynare to solve anything…

Or is this just not going to work? Or have I made a mistake with the timing of the variables? Is it possible to use observation equations(which observe values in the simulation) as variables in the model or would it simply not work in this simple model?

Any help/directions would be greatly appreciated:-)!

P.S I have added the mod file in the appendix, and how the graph is supposed to look
PSV_MA.mod (3.69 KB)

I may be missing something, but when working with

model(linear);
going to third order will not help, because the model is linear and the solution certainty equivalent.