new;
output file=lucas.ot1  reset ; output on;
outwidth 250;
#include c:\gauss\set\ols.set;
#include c:\gauss\set\correlgm.set;

format /rd 7,3;     ts=hsec;

" LUCAS.PG1
 CALIBRATE AND SIMULATE THE RISK PREMIUM IN THE LUCAS MODEL 
 FOR US AND GERMAN DATA.
 TRANSITION PROBABILITIES ESTIMATED BY EVENT COUNTING.

";
@------------------------------------------------------------------@
load X0[97,7]=ch4wrk.prn; @ date, rcons, money@
load spot[97,19]=CH4SPT.prn;    @ Spot exchange rates for 19 countries @

 y0=ln(x0[.,2 3 4]);   @ US, GERMANY, UK  PERCAP REAL CONS@
 m0=ln(x0[.,5 6 7]);   @ US, GERMANY, UK  PERCAP NOMINAL MONEY @
 s0=(-ln(spot[.,7]))~ln(spot[.,2]); @ Germany, UK @


@--------------------- US AND GERMAN DATA ------------------------@
 t0=2; t1=rows(y0);
                             @ US AND GERMAN PER CAPITA MONEY GROWTH RATES @
 dm1=m0[t0:t1,1]-m0[t0-1:t1-1,1];
 dm2=m0[t0:t1,2]-m0[t0-1:t1-1,2];  
                      @ US AND GERMAN GROSS PER CAPITA INCOME GROWTH RATES @
 dy1=1+y0[t0:t1,1]-y0[t0-1:t1-1,1];
 dy2=1+y0[t0:t1,2]-y0[t0-1:t1-1,2];  


t0=1; t1=rows(dy1); sto=zeros(t1,4);
seas=1;
t=t0; do while t le t1;                     @ CREATE SEASONAL DUMMIES @
 if seas eq 1; sto[t,.]=1~0~0~0; endif;
 if seas eq 2; sto[t,.]=0~1~0~0; endif;
 if seas eq 3; sto[t,.]=0~0~1~0; endif;
 if seas eq 4; sto[t,.]=0~0~0~1; endif;
 seas=seas+1; if seas eq 5; seas=1; endif;
t=t+1; endo;
                           @ DESEASONALIZE MONEY BY DUMMY VARIABLE METHOD @
{b,resid1,vcv,tratio,rsq}=ols(dm1,sto,0);
{b,resid2,vcv,tratio,rsq}=ols(dm2,sto,0);

dm1=1+resid1; dm2=1+resid2;     @ GROSS GROWTH RATES @

               @ AVERAGE OF GOOD STATE AND AVERAGE OF BAD STATES @
clear lam1,lam2,lamf1,lamf2,g1,g2,gf1,gf2;
t=1; do while t le t1; @ 1 IS GOOD STATE, 2 IS BAD STATE @
 if dm1[t] gt meanc(dm1); lam1=lam1|dm1[t]; else; lam2=lam2|dm1[t]; endif;
 if dm2[t] gt meanc(dm2); lamf1=lamf1|dm2[t];else;lamf2=lamf2|dm2[t]; endif;
 if dy1[t] gt meanc(dy1); g1=g1|dy1[t]; else; g2=g2|dy1[t]; endif;
 if dy2[t] gt meanc(dy2); gf1=gf1|dy2[t]; else; gf2=gf2|dy2[t]; endif;
t=t+1; endo;

                                           @ GET RID OF FIRST ELT.@
lam1=lam1[2:rows(lam1)];    lam2=lam2[2:rows(lam2)];  
lamf1=lamf1[2:rows(lamf1)]; lamf2=lamf2[2:rows(lamf2)];
g1=g1[2:rows(g1)];          g2=g2[2:rows(g2)];
gf1=gf1[2:rows(gf1)];       gf2=gf2[2:rows(gf2)];

lam1=meanc(lam1);           g1=meanc(g1); 
lam2=meanc(lam2);           g2=meanc(g2); 
lamf1=meanc(lamf1);         gf1=meanc(gf1); 
lamf2=meanc(lamf2);         gf2=meanc(gf2); 


"Average US Consumption Growth good state    :";;g1;
"Average US Consumption Growth bad  state    :";;g2;
"Average German Consumption Growth good state:";;gf1;
"Average German Consumption Growth bad  state:";;gf2;

"Average US Money Growth good state    :";;lam1;
"Average US Money Growth bad  state    :";;lam2;
"Average German Money Growth good state:";;lamf1;
"Average German Money Growth bad  state:";;lamf2;

@---------------------------------------------------------------------@
@::::: ESTIMATE TRANSITION MATRIX ::::::::::::::::@
@---------------------------------------------------------------------@
s=zeros(t1,1);  @ STATE VECTOR @

m_1=meanc(dm1);  m_2=meanc(dm2);
y_1=meanc(dy1);  y_2=meanc(dy2);
ns=zeros(16,1);
t=t0; do while t le t1;
 if dm1[t] gt m_1 and dm2[t] gt m_2 and dy1[t] gt y_1 and dy2[t] gt y_2;  
   s[t]=1; ns[1]=ns[1]+1; endif;
 if dm1[t] gt m_1 and dm2[t] gt m_2 and dy1[t] gt y_1 and dy2[t] le y_2;  
   s[t]=2; ns[2]=ns[2]+1; endif;
 if dm1[t] gt m_1 and dm2[t] gt m_2 and dy1[t] le y_1 and dy2[t] gt y_2;  
   s[t]=3; ns[3]=ns[3]+1; endif;
 if dm1[t] gt m_1 and dm2[t] gt m_2 and dy1[t] le y_1 and dy2[t] le y_2;  
   s[t]=4; ns[4]=ns[4]+1; endif;
 if dm1[t] gt m_1 and dm2[t] le m_2 and dy1[t] gt y_1 and dy2[t] gt y_2;  
   s[t]=5; ns[5]=ns[5]+1; endif;
 if dm1[t] gt m_1 and dm2[t] le m_2 and dy1[t] gt y_1 and dy2[t] le y_2;  
   s[t]=6; ns[6]=ns[6]+1; endif;
 if dm1[t] gt m_1 and dm2[t] le m_2 and dy1[t] le y_1 and dy2[t] gt y_2;  
   s[t]=7; ns[7]=ns[7]+1; endif;
 if dm1[t] gt m_1 and dm2[t] le m_2 and dy1[t] le y_1 and dy2[t] le y_2;  
   s[t]=8; ns[8]=ns[8]+1; endif;

 if dm1[t] le m_1 and dm2[t] gt m_2 and dy1[t] gt y_1 and dy2[t] gt y_2;  
   s[t]=9;  ns[9]=ns[9]+1; endif;
 if dm1[t] le m_1 and dm2[t] gt m_2 and dy1[t] gt y_1 and dy2[t] le y_2;  
   s[t]=10; ns[10]=ns[10]+1; endif;
 if dm1[t] le m_1 and dm2[t] gt m_2 and dy1[t] le y_1 and dy2[t] gt y_2;  
   s[t]=11; ns[11]=ns[11]+1; endif;
 if dm1[t] le m_1 and dm2[t] gt m_2 and dy1[t] le y_1 and dy2[t] le y_2;  
   s[t]=12; ns[12]=ns[12]+1; endif;
 if dm1[t] le m_1 and dm2[t] le m_2 and dy1[t] gt y_1 and dy2[t] gt y_2;  
   s[t]=13; ns[13]=ns[13]+1; endif;
 if dm1[t] le m_1 and dm2[t] le m_2 and dy1[t] gt y_1 and dy2[t] le y_2;  
   s[t]=14; ns[14]=ns[14]+1; endif;
 if dm1[t] le m_1 and dm2[t] le m_2 and dy1[t] le y_1 and dy2[t] gt y_2;  
   s[t]=15; ns[15]=ns[15]+1; endif;
 if dm1[t] le m_1 and dm2[t] le m_2 and dy1[t] le y_1 and dy2[t] le y_2;  
   s[t]=16; ns[16]=ns[16]+1; endif;
t=t+1; endo;


ns[6]=ns[6]-1;   @ I DON'T KNOW WHY, BUT YOU NEED THIS ADJUSTMENT @

p=zeros(16,16);    @ TRANSITION MATRIX @
t=t0+1; do while t le t1;
j=1; do while j le 16;
 if s[t-1] eq j;
  k=1; do while k le 16;
   if s[t] eq k; p[j,k]=p[j,k]+1; endif;
  k=k+1; endo;
 endif;
j=j+1; endo;
t=t+1; endo;

j=1; do while j le 16;
 p[j,.]=p[j,.]/ns[j];
j=j+1; endo;



"Check frequency counts-summing rows of Transition Matrix";;
sumc(p');

format /rd 5,2;
"The transition matrix"; 
p;
format /rd 7,3;
@----------------------------------------------------------------------@
@:::::::::::: INITIAL PROBABILITY VECTOR :::::::::::::::::::::::::::::@
@----------------------------------------------------------------------@
               
 sto=p[1:15,16];
 a=sto~sto~sto~sto~sto~sto~sto~sto~sto~sto~sto~sto~sto~sto~sto;
 a=a-p[1:15,1:15];
 v=inv(eye(15)+a)*sto;
 v=v|(1-sumc(v));    @ THEY ARE EQUALLY LIKELY @

@-------------------------------------------------------------------@
@:::::::::::::: THE LUCAS MODEL ::::::::::::::::::::::::::::::::::::@
@-------------------------------------------------------------------@
theta=0.5; gam=10; beta=0.99;

c=zeros(16,1);  @ THIS IS CORRECT. ERROR FOUND BY MAURICIO BITTENCOURT @
c[1] =(1/lam1)*((g1^(theta))*(gf1^(1-theta)))^(1-gam);
c[2] =(1/lam1)*((g1^(theta))*(gf2^(1-theta)))^(1-gam);
c[3] =(1/lam1)*((g2^(theta))*(gf1^(1-theta)))^(1-gam);
c[4] =(1/lam1)*((g2^(theta))*(gf2^(1-theta)))^(1-gam);
c[5] =(1/lam1)*((g1^(theta))*(gf1^(1-theta)))^(1-gam);
c[6] =(1/lam1)*((g1^(theta))*(gf2^(1-theta)))^(1-gam);
c[7] =(1/lam1)*((g2^(theta))*(gf1^(1-theta)))^(1-gam);
c[8] =(1/lam1)*((g2^(theta))*(gf2^(1-theta)))^(1-gam);
c[9] =(1/lam2)*((g1^(theta))*(gf1^(1-theta)))^(1-gam);
c[10]=(1/lam2)*((g1^(theta))*(gf2^(1-theta)))^(1-gam);
c[11]=(1/lam2)*((g2^(theta))*(gf1^(1-theta)))^(1-gam);
c[12]=(1/lam2)*((g2^(theta))*(gf2^(1-theta)))^(1-gam);
c[13]=(1/lam2)*((g1^(theta))*(gf1^(1-theta)))^(1-gam);
c[14]=(1/lam2)*((g1^(theta))*(gf2^(1-theta)))^(1-gam);
c[15]=(1/lam2)*((g2^(theta))*(gf1^(1-theta)))^(1-gam);
c[16]=(1/lam2)*((g2^(theta))*(gf2^(1-theta)))^(1-gam);


cf=zeros(16,1);
cf[1] =(1/lamf1)*(g1^(theta))*(gf1^(1-theta))^(1-gam);
cf[2] =(1/lamf1)*(g1^(theta))*(gf2^(1-theta))^(1-gam);
cf[3] =(1/lamf1)*(g2^(theta))*(gf1^(1-theta))^(1-gam);
cf[4] =(1/lamf1)*(g2^(theta))*(gf2^(1-theta))^(1-gam);
cf[5] =(1/lamf2)*(g1^(theta))*(gf1^(1-theta))^(1-gam);
cf[6] =(1/lamf2)*(g1^(theta))*(gf2^(1-theta))^(1-gam);
cf[7] =(1/lamf2)*(g2^(theta))*(gf1^(1-theta))^(1-gam);
cf[8] =(1/lamf2)*(g2^(theta))*(gf2^(1-theta))^(1-gam);
cf[9] =(1/lamf1)*(g1^(theta))*(gf1^(1-theta))^(1-gam);
cf[10]=(1/lamf1)*(g1^(theta))*(gf2^(1-theta))^(1-gam);
cf[11]=(1/lamf1)*(g2^(theta))*(gf1^(1-theta))^(1-gam);
cf[12]=(1/lamf1)*(g2^(theta))*(gf2^(1-theta))^(1-gam);
cf[13]=(1/lamf2)*(g1^(theta))*(gf1^(1-theta))^(1-gam);
cf[14]=(1/lamf2)*(g1^(theta))*(gf2^(1-theta))^(1-gam);
cf[15]=(1/lamf2)*(g2^(theta))*(gf1^(1-theta))^(1-gam);
cf[16]=(1/lamf2)*(g2^(theta))*(gf2^(1-theta))^(1-gam);


   @ GROSS HOME CURRENCY DEPRECIATION RATE @
e=zeros(16,1);
e[1]=lam1/lamf1;
e[2]=lam1/lamf1;
e[3]=lam1/lamf1;
e[4]=lam1/lamf1;
e[5]=lam1/lamf2;
e[6]=lam1/lamf2;
e[7]=lam1/lamf2;
e[8]=lam1/lamf2;
e[9] =lam2/lamf1;
e[10]=lam2/lamf1;
e[11]=lam2/lamf1;
e[12]=lam2/lamf1;
e[13]=lam2/lamf2;
e[14]=lam2/lamf2;
e[15]=lam2/lamf2;
e[16]=lam2/lamf2;


                      @ DRAW INITIAL STATE @
seed=418562;
 clear sto;
 u=rndus(1,1,seed);
 if u le v[1]; sto=1; endif;

 if u gt v[1] and u le sumc(v[1:2]) ; sto=2; endif;
 if u gt sumc(v[1:2]) and u  le sumc(v[1:3]); sto=3; endif;
 if u gt sumc(v[1:3]) and u  le sumc(v[1:4]); sto=4; endif;
 if u gt sumc(v[1:4]) and u  le sumc(v[1:5]); sto=5; endif;
 if u gt sumc(v[1:5]) and u  le sumc(v[1:6]); sto=6; endif;
 if u gt sumc(v[1:6]) and u  le sumc(v[1:7]); sto=7; endif;
 if u gt sumc(v[1:7]) and u  le sumc(v[1:8]); sto=8; endif;
 if u gt sumc(v[1:8]) and u  le sumc(v[1:9]); sto=9; endif;
 if u gt sumc(v[1:9]) and u  le sumc(v[1:10]); sto=10; endif;
 if u gt sumc(v[1:10]) and u le sumc(v[1:11]); sto=11; endif;
 if u gt sumc(v[1:11]) and u le sumc(v[1:12]); sto=12; endif;
 if u gt sumc(v[1:12]) and u le sumc(v[1:13]); sto=13; endif;
 if u gt sumc(v[1:13]) and u le sumc(v[1:14]); sto=14; endif;
 if u gt sumc(v[1:14]) and u le sumc(v[1:15]); sto=15; endif;
 if u gt sumc(v[1:15]) and u le 1; sto=16; endif;

state=sto;   @ INITIAL STATE @

@-------------------------------------------------------------@
@::::::::::: COMPUTE STATE-CONTINGENT RISK PREMIA ::::::::::::@
@-------------------------------------------------------------@
rp=zeros(16,1);   @ 16 POSSIBLE STATES.  @

k=1; do while k le 16;
 rp[k]=p[k,.]*e-(p[k,.]*cf)/(p[k,.]*c);
k=k+1; endo;

"State contingent risk premia in percent";
seqa(1,1,16)~rp*100;

@-------------------------------------------------------------@
@:::::::::::::: Generate a realization :::::::::::::::::::::::@
@-------------------------------------------------------------@

t1=500;

t=2; do while t le t1;  @ GENERATE A SEQUENCE OF STATES @
 k=state[rows(state)];  @ CURRENT STATE @
 clear sto;
 u=rndus(1,1,seed);
  if u le p[k,1]; sto=1; endif;
  if u gt p[k,1] and u le sumc(p[k,1:2]'); sto=2; endif;
  if u gt sumc(p[k,1:2]') and u  le sumc(p[k,1:3]'); sto=3; endif;
  if u gt sumc(p[k,1:3]') and u  le sumc(p[k,1:4]'); sto=4; endif;
  if u gt sumc(p[k,1:4]') and u  le sumc(p[k,1:5]'); sto=5; endif;
  if u gt sumc(p[k,1:5]') and u  le sumc(p[k,1:6]'); sto=6; endif;
  if u gt sumc(p[k,1:6]') and u  le sumc(p[k,1:7]'); sto=7; endif;
  if u gt sumc(p[k,1:7]') and u  le sumc(p[k,1:8]'); sto=8; endif;
  if u gt sumc(p[k,1:8]') and u  le sumc(p[k,1:9]'); sto=9; endif;
  if u gt sumc(p[k,1:9]') and u  le sumc(p[k,1:10]'); sto=10; endif;
  if u gt sumc(p[k,1:10]') and u le sumc(p[k,1:11]'); sto=11; endif;
  if u gt sumc(p[k,1:11]') and u le sumc(p[k,1:12]'); sto=12; endif;
  if u gt sumc(p[k,1:12]') and u le sumc(p[k,1:13]'); sto=13; endif;
  if u gt sumc(p[k,1:13]') and u le sumc(p[k,1:14]'); sto=14; endif;
  if u gt sumc(p[k,1:14]') and u le sumc(p[k,1:15]'); sto=15; endif;
  if u gt sumc(p[k,1:15]') and u le 1;                sto=16; endif;

state=state|sto; clear sto;
t=t+1; endo;

slev=zeros(t1,1); @ Exchange Rate levels @
sg=zeros(t1,1);  @ GROSS GROWTH RATE IN SPOT EXCHANGE RATE @
fp=zeros(t1,1);  @ FORWARD PREMIUM @
erp=zeros(t1,1); @ REALIZED  RISK PRMEMIUM @
t=1; do while t le t1;
 sg[t]=e[state[t]];
/*------------------------------------------------------------
 bf=beta*p[state[t],.]*c;      @<=== Xiaodai Xiu found this mistake @
 b=beta*p[state[t],.]*cf;
------------------------------------------------------*/
 bf=beta*p[state[t],.]*cf;
 b=beta*p[state[t],.]*c;

 fp[t]=bf/b;
 erp[t]=rp[state[t]];
t=t+1; endo;

slev[1]=1;
t=2; do while t le t1;
 slev[t]=slev[t-1]*sg[t];
t=t+1; endo;

@--------------------------------------------------------------@

library pgraph;
graphset;
      fonts("simplex complex microb simgrma");   @ LOAD FONTS @
@     _pline = { 1 6 55 0 1995 0 1 2 1,    @
@                1 6 55 1 1995 1 1 2 1} ;  @
     _pltype = { 6  }; /* line type for the main curves */
     _psymsiz = {3};      @ SYMBOL SIZE @
     _plwidth = {5};      @ LINE THICKNESS  @
     _pdate="";
  xtics(73,97,2,1);
 _plctrl = {0 1};  @ First line solid. Second line symbol every 1 pt @
 _ptitlht = {0.2};   @ TITLE SIZE @
 _pstype = {9};
 _plegstr = "";
@----------------------------------------------------------@
caltme=seqa(73,0.25,rows(fp));
begwind; window(1,2,0);
title("\202A"); @ Tomorrow's gross depreciation (S_{t+1}/S_t)@
                @ and today's forward premium F_t/S_t (boxes) @
xy(caltme[2:t1],sg[2:t1]~fp[1:t1-1]);

 nextwind;
title("\202B");
           @ Tomorrow's realized speculative profit @
           @ and today's expected excess return (risk premiuum in boxes) @
xy(caltme[2:t1],(sg[2:t1]-fp[1:t1-1])~erp[1:t1-1]);
endwind;

format /rd 15,4; 

@xy(caltme[2:t1],slev[2:t1]);@

" Write out observations for plotting in Excel";
@format /re 14,6;@
caltme[2:97]~sg[2:97]~fp[1:96]~caltme[2:97]~(sg[2:97]-fp[1:97-1])~erp[1:97-1]

@----------------------------------------------------------@
"Regression of future depreciation on forward premium ";
{b,resid,vcv,tratio,rsq}=ols(ln(sg[2:t1]),ln(fp[1:t1-1]),1);
"Coeff  t-ratio";
b~tratio;
correlgm(resid,10);

"Regression of S'/S on F/S (gross depreciation on F/S) ";
{b,resid,vcv,tratio,rsq}=ols(sg[2:t1],fp[1:t1-1],1);
"Coeff  t-ratio";
b~tratio;

@----------------------------------------------------------@

"Volatility of risk premium            ";; stdc(rp);
"Volatility of realized profit (F-S')/S";; stdc(e-c./cf);
"Volatility of forward premium (F/S)   ";; stdc(fp);
"Volatility of depreciation (S'/S)     ";; stdc(sg);

"Autocorrelation of risk premium       ";; correlgm(rp,2)';
"Autocorrelation of realized profit    ";; correlgm(e-c./cf,2)';
"Autocorrelation of forward premium    ";; correlgm(fp,2)';
"Autocorrelation of depreciation       ";; correlgm(sg,2)';

output off;