stepan
/
dynare
1
0
Fork 0
Code Issues Pull Requests Packages Projects Releases Wiki Activity

1st cut of KF refactoring and simplification 

git-svn-id: https://www.dynare.org/svn/dynare/trunk@2750 ac1d8469-bf42-47a9-8791-bf33cf982152
time-shift
george 2009-06-10 06:40:07 +00:00
parent 36ecf9d6bd
commit 1245d60326
5 changed files with 1836 additions and 1557 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

189
mex/sources/kalman/cc/kalman.cpp
View File

@ -468,6 +468,52 @@ void SmootherResults::exportV(GeneralMatrix&out)const
}
}
BasicKalmanTask::BasicKalmanTask(const GeneralMatrix&d,const GeneralMatrix&ZZ,
const GeneralMatrix&HH,const GeneralMatrix&TT,
const GeneralMatrix&RR,const GeneralMatrix&QQ,
const StateInit&init_state)
: // ssf(Z,H,T,R,Q),
data(d), Zt(*(new ConstGeneralMatrix(ZZ))),
Ht(*(new ConstGeneralMatrix(HH))),
Tt(*(new ConstGeneralMatrix(TT))),
Rt(*(new ConstGeneralMatrix(RR))),
Qt(*(new ConstGeneralMatrix(QQ))),
init(init_state)
{
TS_RAISE_IF(d.numRows()!=Zt.numRows(),
"Data not compatible with BasicKalmanTask constructor");
// TS_RAISE_IF(ssf.m!=init.getM(),
// "State initialization not compatible with SSF in KalmanTask constructor");
}
BasicKalmanTask::BasicKalmanTask(const GeneralMatrix&d,const ConstGeneralMatrix&ZZ,
const ConstGeneralMatrix&HH,const ConstGeneralMatrix&TT,
const ConstGeneralMatrix&RR,const ConstGeneralMatrix&QQ,
const StateInit&init_state)
: // ssf(Z,H,T,R,Q),
data(d), Zt(ZZ), Ht(HH), Tt(TT), Rt(RR), Qt(QQ),init(init_state)
{
TS_RAISE_IF(d.numRows()!=Zt.numRows(),
"Data not compatible with BasicKalmanTask constructor");
// TS_RAISE_IF(ssf.m!=init.getM(),
// "State initialization not compatible with SSF in KalmanTask constructor");
}
BasicKalmanTask::~BasicKalmanTask()
{
if (&Zt);
delete &Zt;
if (&Ht);
delete &Ht;
if (&Tt);
delete &Tt;
if (&Rt);
delete &Rt;
if (&Qt);
delete &Qt;
}
KalmanTask::KalmanTask(const GeneralMatrix&d,const GeneralMatrix&Z,
const GeneralMatrix&H,const GeneralMatrix&T,
const GeneralMatrix&R,const GeneralMatrix&Q,
@ -556,6 +602,13 @@ KalmanTask::filter(int&per,int&d, int start, std::vector<double>* vll)const
return 0.0;
}
double
BasicKalmanTask::filter(int&per,int&d, int start, std::vector<double>* vll)const
{
d= 0;
per= vll->size() ;
return filterNonDiffuse(init.getA(),init.getPstar(), start, vll);
}
/*****************
This is public interface that runs a filter followed by a smoother.
In addition to things returned by |KalmanTask::filter|, it fills also
@ -594,6 +647,139 @@ KalmanTask::filter_and_smooth(SmootherResults&sres,
return 0.0;
}
/*****************
This runs a Basic non-diffuse filter with the given $t$, $a_t$ and
$P_t$. It fills the |FilterResults|.
First we check that the passed $P_t$ is positive definite by checking
that it has not negative diagonal and is symmetric diagonally
dominant. This is not equivalent to positive definitness but it excludes
``much'' of indefinite matrices. This check is important since this
routine is called from a diffuse filter and it is possible due to a
wrong guess/calculation of a number of diffuse periods that numerical
instability and roundoff errors make the matrix $P_*$ broken.
Then we cycle until the end of period and perform a classical Kalman
filter operations.
*****************/
double
BasicKalmanTask::filterNonDiffuse(const Vector&a,const GeneralMatrix&P,
int start, std::vector<double>* vll)const
{
double loglik=0;
Vector at(a);
GeneralMatrix Pt(P);
if(TSUtils::hasNegativeDiagonal(Pt)||!TSUtils::isSymDiagDominant(Pt))
return 0.0;
/*
ConstGeneralMatrix Zt(Z);
ConstGeneralMatrix Ht(H);
ConstGeneralMatrix Tt(T);
ConstGeneralMatrix Qt(Q);
ConstGeneralMatrix Rt(R);
*/
bool isTunit=0;// Tt->isUnit();
bool isQzero= Qt.isZero();
bool isRzero= Rt.isZero();
int t= 1;
int nonsteady=1;
for(;t<=data.numCols()&&nonsteady;++t)
{
ConstVector yt(data,t-1);
/*****************
This calculates $$v_t = y_t - Z_t*a_t.$$
*****************/
Vector vt(yt);
Zt.multsVec(vt,at);
/*****************
This calculates $$F_t = Z_tP_tZ_t^T+H_t.$$
*****************/
GeneralMatrix Mt(Pt,Zt,"trans");
GeneralMatrix Ft(Ht);
Ft.multAndAdd(Zt,ConstGeneralMatrix(Mt));
PLUFact Ftinv(Ft);
if(!Ftinv.isRegular())
return 0.0;
/*****************
This calculates $$K_t = T_tP_tZ_t^TF_t^{-1}.$$
*****************/
GeneralMatrix Kt(Tt,Mt);
Ftinv.multInvRight(Kt);
/*****************
This calculates $$L_t = T_t-K_tZ_t.$$
*****************/
GeneralMatrix Lt(Tt);
Lt.multAndAdd(ConstGeneralMatrix(Kt),Zt,-1.0);
/*****************
Here we calc likelihood and store results.
*****************/
double ll= calcStepLogLik(Ftinv,vt);
// fres.set(t,Ftinv,vt,Lt,at,Pt,ll);
(*vll)[t-1]=ll;
if (t>start) loglik+=ll;
if(t<data.numCols())
{
/*****************
This calculates $$a_{t+1} = T_ta_t + K_tv_t.$$
*****************/
if(!isTunit)
{
Vector atsave((const Vector&)at);
Tt.multVec(0.0,at,1.0,atsave);
}
Kt.multVec(1.0,at,1.0,ConstVector(vt));
/*****************
This calculates $$P_{t+1} = T_tP_tL_t^T + R_tQ_tR_t^T.$$
*****************/
GeneralMatrix PtLttrans(Pt,Lt,"trans");
if(!isTunit)
{
Pt.zeros();
Pt.multAndAdd(Tt,ConstGeneralMatrix(PtLttrans));
}
else
{
Pt= (const GeneralMatrix&)PtLttrans;
}
if(!isRzero&&!isQzero)
{
GeneralMatrix QtRttrans(Qt,Rt,"trans");
Pt.multAndAdd(Rt,ConstGeneralMatrix(QtRttrans));
}
}
}
// for(;t<=data.numCols();t++)
// {
// ConstVector yt(data,t-1);
// }
return loglik;
}
double
BasicKalmanTask::calcStepLogLik(const PLUFact&Finv,const Vector&v)
{
int p= Finv.numRows();
Vector Finvv(v);
Finv.multInvLeft(Finvv);
double vFinvv= v.dot(Finvv);
return-0.5*(p*log(2*M_PI)+Finv.getLogDeterminant()+vFinvv);
}
/*****************
This runs a non-diffuse filter with the given $t$, $a_t$ and
$P_t$. It fills the |FilterResults|.
@ -609,7 +795,6 @@ instability and roundoff errors make the matrix $P_*$ broken.
Then we cycle until the end of period and perform a classical Kalman
filter operations.
*****************/
void
KalmanTask::filterNonDiffuse(const Vector&a,const GeneralMatrix&P,
int first,FilterResults&fres)const
@ -618,6 +803,7 @@ KalmanTask::filterNonDiffuse(const Vector&a,const GeneralMatrix&P,
GeneralMatrix Pt(P);
if(TSUtils::hasNegativeDiagonal(Pt)||!TSUtils::isSymDiagDominant(Pt))
return;
for(int t= first;t<=data.numCols();t++)
{
ConstVector yt(data,t-1);
@ -1316,6 +1502,7 @@ This evaluates a step loglikelihood as
v_t^TF_t^{-1}v_t\right]$$
This is a static method.
*****************/
double
KalmanTask::calcStepLogLik(const PLUFact&Finv,const Vector&v)
{

44
mex/sources/kalman/cc/kalman.h
View File

@ -154,6 +154,50 @@ class SmootherResults{
};
class BasicKalmanTask{
// friend class KalmanUniTask;
// SSForm ssf;
const GeneralMatrix &data;
const ConstGeneralMatrix &Zt;
const ConstGeneralMatrix &Ht;
const ConstGeneralMatrix &Tt;
const ConstGeneralMatrix &Rt;
const ConstGeneralMatrix &Qt;
const StateInit&init;
public:
BasicKalmanTask(const GeneralMatrix&d,const GeneralMatrix&ZZ,
const GeneralMatrix&HH,const GeneralMatrix&TT,
const GeneralMatrix&RR,const GeneralMatrix&QQ,
const StateInit&init_state);
// BasicKalmanTask(const GeneralMatrix&d,const TMatrix&Z,
// const TMatrix&H,const TMatrix&T,
// const TMatrix&R,const TMatrix&Q,
// const StateInit&init_state);
BasicKalmanTask(const GeneralMatrix&d,const ConstGeneralMatrix&ZZ,
const ConstGeneralMatrix&HH,const ConstGeneralMatrix&TT,
const ConstGeneralMatrix&RR,const ConstGeneralMatrix&QQ,
const StateInit&init_state);
virtual ~BasicKalmanTask();
// double filter(int&per,int&d)const;
// double filter(int&per,int&d, int start, std::vector<double>* vll)const;
double filter(int&per,int&d,int start, std::vector<double>* vll)const;
// double filter_and_smooth(SmootherResults&sres,int&per,int&d)const;
protected:
double filterNonDiffuse(const Vector&a,const GeneralMatrix&Pstar,
int start, std::vector<double>* vll) const; //int first,FilterResults&fres)const;
// void filterDiffuse(const Vector&a,const GeneralMatrix&Pstar,
// const GeneralMatrix&Pinf,int first,
// DiffuseFilterResults&fres)const;
// void smootherNonDiffuse(const FilterResults&fres,SmootherResults&sres)const;
// void smootherDiffuse(const DiffuseFilterResults&fres,SmootherResults&sres)const;
// void smootherNonDiffuseStep(int t,const FilterResults&fres,
// Vector&rt,GeneralMatrix&Nt,
// Vector&alphat,GeneralMatrix&Vt,
// Vector&etat)const;
static double calcStepLogLik(const PLUFact&Finv,const Vector&v);
};
class KalmanUniTask;
class KalmanTask{
friend class KalmanUniTask;

28
mex/sources/kalman/matlab/Makefile
View File

@ -1,16 +1,16 @@
# $Id: Makefile 531 2005-11-30 13:49:48Z kamenik $
# Copyright 2005, Ondra Kamenik
DEBUG = yes
#DEBUG = yes
#LD_LIBS := -llapack -lcblas -lf77blas -latlas -lg2c
CC_FLAGS := -DMATLAB -DWINDOWS -DNO_BLAS_H -DNO_LAPACK_H \
-Wall -I../../gensylv/cc -I../cc \
-Wall -I../sylv/cc -I../cc \
-Ic:/"Program Files"/MATLAB_SV71/extern/include
ifeq ($(DEBUG),yes)
# CC_FLAGS := -DDEBUG $(CC_FLAGS) -g
CC_FLAGS := -DTIMING_LOOP -DDEBUG $(CC_FLAGS) -g -pg #-Wl,-pg
CC_FLAGS := -DTIMING_LOOP -DDEBUG $(CC_FLAGS) -g #-pg #-Wl,-pg
KALMANLIB := kalmanlib_dbg.a
else
# CC_FLAGS := $(CC_FLAGS) -O2
@ -34,9 +34,9 @@ endif
# end add
matrix_interface := GeneralMatrix Vector SylvException
matobjs := $(patsubst %, ../../gensylv/cc/%.o, $(matrix_interface))
mathsource := $(patsubst %, ../../gensylv/cc/%.h, $(matrix_interface))
matcppsource := $(patsubst %, ../../gensylv/cc/%.cpp, $(matrix_interface))
matobjs := $(patsubst %, ../sylv/cc/%.o, $(matrix_interface))
mathsource := $(patsubst %, ../sylv/cc/%.h, $(matrix_interface))
matcppsource := $(patsubst %, ../sylv/cc/%.cpp, $(matrix_interface))
# cppsource := $(patsubst %.cweb,%.cpp,$(cwebsource))
kalmancppsource := $(wildcard ../cc/*.cpp)
kalmanhsource := $(wildcard ../cc/*.h)
@ -62,27 +62,27 @@ dummy.ch:
c++ $(CC_FLAGS) -c $*.cpp
kalmanlib.a: $(objects) # $(matobjs) #$(kalmanobjects)
ar cr kalmanlib.a $(kalmanobjects) $(matobjs) # $(objects)
ranlib kalmanlib.a
$(KALMANLIB): $(objects) # $(matobjs) #$(kalmanobjects)
ar cr $(KALMANLIB) $(kalmanobjects) $(matobjs) # $(objects)
ranlib $(KALMANLIB)
kalman_smoother_dll.dll: kalman_smoother.o kalmanlib.a #$(hsource) $(cppsource)
kalman_smoother_dll.dll: kalman_smoother.o $(KALMANLIB) #$(hsource) $(cppsource)
gcc -shared $(CC_FLAGS) -o kalman_smoother_dll.dll kalman_smoother.o \
kalmanlib.a $(LD_LIBS)
minv.dll: minv.o kalmanlib.a # $(hsource) $(cppsource)
minv.dll: minv.o $(KALMANLIB) # $(hsource) $(cppsource)
gcc -shared $(CC_FLAGS) -o minv.dll minv.o \
kalmanlib.a $(LD_LIBS)
kalman_filter_dll.dll: kalman_filters.o kalmanlib.a # $(hsource) $(cppsource)
kalman_filter_dll.dll: kalman_filters.o $(KALMANLIB) # $(hsource) $(cppsource)
gcc -shared $(CC_FLAGS) -o kalman_filter_dll.dll kalman_filters.o \
$(KALMANLIB) $(LD_LIBS)
kalman_filters_testx.exe: kalman_filters_testx.o kalmanlib.a # $(hsource) $(cppsource)
kalman_filters_testx.exe: kalman_filters_testx.o $(KALMANLIB) # $(hsource) $(cppsource)
gcc $(CC_FLAGS) -pg -o kalman_filters_testx.exe kalman_filters_testx.o \
$(KALMANLIB) $(LD_LIBS)
all: $(objects) kalmanlib.a kalman_smoother_dll.dll kalman_filter_dll.dll # $(cppsource) $(hsource) $(kalmanhsource) $(kalmancppsource)
all: $(objects) $(KALMANLIB) kalman_smoother_dll.dll kalman_filter_dll.dll # $(cppsource) $(hsource) $(kalmanhsource) $(kalmancppsource)
kalman_filter_loop.o: kalman_filters.cpp
c++ -DTIMING_LOOP $(CC_FLAGS) -o kalman_filter_loop.o kalman_filters.cpp

26
mex/sources/kalman/matlab/kalman_filters.cpp
View File

@ -115,7 +115,10 @@ extern "C" {
// create state init
StateInit* init = NULL;
std::vector<double>* vll=new std::vector<double> (nper);
if (diffuse) {
if (diffuse||uni)
{
if (diffuse)
{
GeneralMatrix Pinf(mxGetPr(prhs[9]), mxGetM(prhs[9]), mxGetN(prhs[9]));
init = new StateInit(P, Pinf, a.getData());
}
@ -141,9 +144,28 @@ extern "C" {
loglik = kt.filter(per, d, (start-1), vll);
#ifdef TIMING_LOOP
}
//mexPrintf("kalman_filter: finished 10,000 loops");
mexPrintf("kalman_filter: finished 1000 loops");
#endif
}
}
else // basic Kalman
{
init = new StateInit(P, a.getData());
BasicKalmanTask bkt(Y, Z, H, T, R, Q, *init);
#ifdef TIMING_LOOP
for (int tt=0;tt<1000;++tt)
{
#endif
loglik = bkt.filter( per, d, (start-1), vll);
#ifdef DEBUG
mexPrintf("Basickalman_filter: loglik=%d \n", loglik);
#endif
#ifdef TIMING_LOOP
}
mexPrintf("Basickalman_filter: finished 1000 loops");
#endif
}
// destroy init
delete init;

28
mex/sources/kalman/matlab/kalman_filters_testx.cpp
View File

@ -219,7 +219,11 @@ int main(int argc, char* argv[])
// create state init
StateInit* init = NULL;
std::vector<double>* vll=new std::vector<double> (nper);
if (diffuse) {
if (diffuse||uni)
{
if (diffuse)
{
GeneralMatrix Pinf(P.numRows(),P.numCols());
Pinf.zeros();
init = new StateInit(P, Pinf, a.getData());
@ -250,6 +254,28 @@ int main(int argc, char* argv[])
}
mexPrintf("kalman_filter: finished 10,000 loops");
#endif
}
else // basic Kalman
{
init = new StateInit(P, a.getData());
BasicKalmanTask bkt(Y, Z, H, T, R, Q, *init);
#ifdef TIMING_LOOP
for (int tt=0;tt<10000;++tt)
{
#endif
loglik = bkt.filter( per, d, (start-1), vll);
#ifdef DEBUG
mexPrintf("Basickalman_filter: loglik=%d \n", loglik);
#endif
#ifdef TIMING_LOOP
}
mexPrintf("Basickalman_filter: finished 10,000 loops");
#endif
}
// destroy init
delete init;
mexPrintf("logLik = %f \n", loglik);