[Tests] Fix RuntimeWarning in lqcontrol.py

[oyamad]

/home/travis/build/QuantEcon/QuantEcon.py/quantecon/lqcontrol.py:244: RuntimeWarning: invalid value encountered in double_scalars
  d = self.beta * np.trace(dot(P, dot(C, C.T))) / (1 - self.beta)

[duncanhobbs]

I think this is a divide by zero error. self.beta (the discount parameter) is set to 1 by default.

Does it make sense to shift the discount parameter default to something like 0.99?

[oyamad]

Thanks @duncanhobbs

When T==None (i.e., the model is one of infinite horizon), then beta=1 should not be accepted.

It should be checked here:

QuantEcon.py/quantecon/lqcontrol.py

Lines 137 to 146 in d9c1e8a

	if T:
	# == Model is finite horizon == #
	self.T = T
	self.Rf = np.asarray(Rf, dtype='float')
	self.P = self.Rf
	self.d = 0
	else:
	self.P = None
	self.d = None
	self.T = None

The pair of default values (beta=1, T=None) is not a good choice:

QuantEcon.py/quantecon/lqcontrol.py

Line 115 in d9c1e8a

def __init__(self, Q, R, A, B, C=None, N=None, beta=1, T=None, Rf=None):

@jstac ?

[jstac]

I'm a bit confused. I think beta=1 is quite common in infinite horizon LQ models, since costs can be driven to zero.

See, e.g.,

https://en.wikipedia.org/wiki/Linear%E2%80%93quadratic_regulator#Infinite-horizon,_discrete-time_LQR

Yet it clearly causes an issue here. I have to think about this when I have a bit of time.

[jstac]

Perhaps we can't have beta=1 in infinite horizons when there's a stochastic component to the update rule for the state, since it keeps pushing the state away from zero and hence pushing up the cost.

So first best is that we allow beta=1 when C=0 and force it to be < 1 otherwise. Agree?

[oyamad]

@jstac Thanks, but what if C != 0 and T != None (i.e., finite horizon)? (nevermind)

And here (in the method stationary_values which is for infinite horizon)

QuantEcon.py/quantecon/lqcontrol.py

Line 244 in d9c1e8a

d = self.beta * np.trace(dot(P, dot(C, C.T))) / (1 - self.beta)

it should be something like:

if self.beta == 1:
    d = 0.
else:
    d = self.beta * np.trace(dot(P, dot(C, C.T))) / (1 - self.beta)

[oyamad]

Fixed by #509.