/*
**  Durat0
**
**  Purpose:
**    Generate data from duration model
**
**  Date:
**    10/8/09
**
**  Author:
**    Charles Bos
*/
#include <oxstd.h>
#include <oxprob.h>
#include <oxdraw.h>

/*
**  GenerateData(const avY, const amX, const vBeta, const dAlpha, const iN)
**
**  Purpose:
**    Generate both exogenous variables and duration data
**
**  Inputs:
**    vBeta     iK x 1 vector of regression-type parameters
**    dAlpha    double, duration dependence
**    iN        integer, number of observations
**
**  Outputs:
**    avY       address of iN x 1 vector of durations
**    amX       address of iN x iK matrix with exogenous variables,
**              containing a constant and remaining N(0,1) columns
**
**  Return value:
**    none
*/
GenerateData(const avY, const amX, const vBeta, const dAlpha, const iN)
{
  decl vL, vE, iK;
  
  iK= sizerc(vBeta);
  amX[0]= 1~rann(iN, iK-1);
  vL= exp(amX[0]*vBeta);
  
  vE= ranexp(iN, 1, 1.0);
  avY[0]= (vE.^(1/dAlpha)) ./ vL;
}

/*
**  Output(const vY, const mX, const vBeta, const dAlpha)
**
**  Purpose:
**    Present some output on the data 
**
**  Inputs:
**    vY        iN x 1 vector of durations
**    mX        iN x iK matrix with exogenous variables,
**              containing a constant and remaining N(0,1) columns
**    vBeta     iK x 1 vector of regression-type parameters
**    dAlpha    double, duration dependence
**
**  Return value:
**    none
*/
Output(const vY, const mX, const vBeta, const dAlpha)
{
  decl mR;
  
  mR= vY~mX;
  println ("Generated ", sizerc(vY), 
           " observations from Weibull with beta=", vBeta', 
           "and alpha= ", dAlpha);
  
  println ("Mean and standard deviations: ", 
           "%c", {"y", "x0", "x1"},
           "%r", {"m", "s"},
           meanc(mR)|sqrt(varc(mR)));
           
  DrawDensity(0, vY', {"y"}, 1, 0, 0);
  DrawMatrix(1, vY', "y", 1, 1);
  ShowDrawWindow();           
}


main()
{
  decl iN, vY, mX, vL, vBeta, dAlpha, sOut;
  
  // Magic numbers 
  iN= 1000;
  vBeta= <1; 1>;
  dAlpha= 1.5;
  sOut= "data/genrdur.fmt";
  
  GenerateData(&vY, &mX, vBeta, dAlpha, iN);
  Output(vY, mX, vBeta, dAlpha);
  
  savemat(sOut, vY~mX);  

}