Question about replication of Filardo(1994)

[fan]

Dear Tom,

I am trying to replicate Filardo(1994) by using MLE method. I am keeping getting the following error message

## MAT15. Subscripts Too Large or Non-Positive
The Error Occurred At Location 216, Line 21 of %MSVARPROB

I am not sure how to solve the problem and I am seeking for your help.

Here is my code

Code: Select all

cal(m) 1948
open data filardo.dat
data(format=prn,org=cols) 1948:1 1992:8 year month ip cli dfi xli sp fed tspread
*
set ipgr = 100.0*log(ip/ip{1})
*
* The graph in the paper appears to be standardized (probably by dividing through
* by 2.0 times the sample standard deviation).
*
graph(footer="Figure 1. Log Growth Rate of Monthly IP")
# ipgr
*
stats(noprint) ipgr * 1959:12
compute stdearly=sqrt(%variance)
stats(noprint) ipgr 1960:1 *
compute stdlate=sqrt(%variance)
*
set ipgradjust = %if(t<=1959:12,ipgr*stdlate/stdearly,ipgr)
*
set g = ipgradjust
*
* Transform the information variables
*
set cligr   = 100.0*log(cli/cli{1})
set spgr    = 100.0*log(sp/sp{1})
set xlidiff = xli-xli{1}
set ffdiff  = fed-fed{1}
set tspdiff = tspread-tspread{1}
set dfismooth = dfi+2*dfi{1}+2*dfi{2}+dfi{3}
*
* And center them to zero means
*
dofor s = cligr spgr xlidiff ffdiff tspdiff dfismooth
   diff(center) s
end dofor s
*******************************************************
*
* Common setup
*
compute nlags=4
*
source msvarsetup.src
compute gstart=1948:6,gend=1992:8
@MSVARSetup(lags=nlags)
# g
*********************************************************************
*
*
nonlin(parmset=common) mu phi sigma

*********************************************************************
*
* Time varying probabilities
*
dec equation p1eq p2eq
dec vector   v1   v2
nonlin(parmset=tv) v1 v2
*********************************************************************
*
* This overrides the standard fixed definition
*
function  %MSVARPmat  time
type rect %MSVARPmat
type int  time
*
local integer i j
local rect    pexpand
local real    z
*
compute p(1,1)=%(z=exp(%dot(%eqnxvector(p1eq,time),v1)),z/(1+z))
compute p(1,2)=%(z=exp(%dot(%eqnxvector(p2eq,time),v2)),1/(1+z))
compute %MSVARInitTransition()
*
if nexpand==nstates {
   dim pexpand(nstates,nstates-1)
   ewise pexpand(i,j)=%if(i==nstates,-1.0,i==j)
   compute %MSVARPmat=pexpand*p
   ewise %MSVARPmat(i,j)=%if(i==nstates,1.0+%MSVARPmat(i,j),%MSVARPmat(i,j))
}
else {
   dim %MSVARPmat(nexpand,nexpand)
   ewise %MSVARPmat(i,j)=MSVARTransProbs(MSVARTransLookup(i,j))
}
end
*********************************************************************
*
* Initialize the probabilities using only the intercepts
*
function TVPInit
type vector TVPInit
*
local rect a
local integer i j
*
compute p=||%logistic(v1(1),1.0),1-%logistic(v2(1),1.0)||
compute TVPInit=%MSVARInit()
end
*********************************************************************
frml logltvtp = %MSVARPMat(t),log(%MSVarProb(t))
*********************************************************************
*
* Define the logistic index for the transitions
*
equation p1eq *
# constant cligr{1}
equation p2eq *
# constant cligr{1}
*
* Standard guess values aiming to make state 1 the low growth state.
*
compute v1=log(.7/.3)~~0.0
compute v2=log(.95/.05)~~0.0
compute mu(1)=-1.7,mu(2)=.3

boxjenk(ar=nlags) g
compute phi=||%beta(1),%beta(2),%beta(3),%beta(4)||
compute sigma=%seesq
*
maximize(parmset=tv+common,start=(pstar=TVPInit()),$
 method=bhhh,iters=20) logltvtp gstart gend

Thank you in advance

[TomDoan]

Add an

@MSFilterInit

before the MAXIMIZE. It makes sure everything is set up properly.

[fan]

quote="TomDoan"]Add an

@MSFilterInit

before the MAXIMIZE. It makes sure everything is set up properly.[/quote]

Thank you for the quick reply. Since @MSFliterInit can't be used together with msvarsetup, I changed my code to mssetup. However, after setting up the ms procedure, the new error message says [quote] ## I2. Expected Instruction Here
>>>>#<<<<[quote] I could not figure out the mistake I made, so I am seeking for your help.

Here is my modified code

Code: Select all

cal(m) 1948 a
open data filardo.dat
data(format=prn,org=cols) 1948:1 1992:8 year month ip cli dfi xli sp fed tspread
*
set ipgr = 100.0*log(ip/ip{1})
*
* The graph in the paper appears to be standardized (probably by dividing through
* by 2.0 times the sample standard deviation).
*
graph(footer="Figure 1. Log Growth Rate of Monthly IP")
# ipgr
*
stats(noprint) ipgr * 1959:12
compute stdearly=sqrt(%variance)
stats(noprint) ipgr 1960:1 *
compute stdlate=sqrt(%variance)
*
set ipgradjust = %if(t<=1959:12,ipgr*stdlate/stdearly,ipgr)
*
set g = ipgradjust
*
* Transform the information variables
*
set cligr   = 100.0*log(cli/cli{1})
set spgr    = 100.0*log(sp/sp{1})
set xlidiff = xli-xli{1}
set ffdiff  = fed-fed{1}
set tspdiff = tspread-tspread{1}
set dfismooth = dfi+2*dfi{1}+2*dfi{2}+dfi{3}
*
* And center them to zero means
*
dofor s = cligr spgr xlidiff ffdiff tspdiff dfismooth
   diff(center) s
end dofor s
*******************************************************
*
* Common setup
*
compute nlags=4
*
source msvarsetup.src
compute gstart=1948:6,gend=1992:8
@MSVARSetup(lags=nlags)
# g
*********************************************************************
*
*
nonlin(parmset=common) mu phi sigma

*********************************************************************
*
* Time varying probabilities
*
dec equation p1eq p2eq
dec vector   v1   v2
nonlin(parmset=tv) v1 v2
*********************************************************************
*
* This overrides the standard fixed definition
*
function  %MSVARPmat  time
type rect %MSVARPmat
type int  time
*
local integer i j
local rect    pexpand
local real    z
*
compute p(1,1)=%(z=exp(%dot(%eqnxvector(p1eq,time),v1)),z/(1+z))
compute p(1,2)=%(z=exp(%dot(%eqnxvector(p2eq,time),v2)),1/(1+z))
compute %MSVARInitTransition()
*
if nexpand==nstates {
   dim pexpand(nstates,nstates-1)
   ewise pexpand(i,j)=%if(i==nstates,-1.0,i==j)
   compute %MSVARPmat=pexpand*p
   ewise %MSVARPmat(i,j)=%if(i==nstates,1.0+%MSVARPmat(i,j),%MSVARPmat(i,j))
}
else {
   dim %MSVARPmat(nexpand,nexpand)
   ewise %MSVARPmat(i,j)=MSVARTransProbs(MSVARTransLookup(i,j))
}
end
*********************************************************************
*
* Initialize the probabilities using only the intercepts
*
function TVPInit
type vector TVPInit
*
local rect a
local integer i j
*
compute p=||%logistic(v1(1),1.0),1-%logistic(v2(1),1.0)||
compute TVPInit=%MSVARInit()
end
*********************************************************************
frml logltvtp = %MSVARPMat(t),log(%MSVarProb(t))
*********************************************************************
*
* Define the logistic index for the transitions
*
equation p1eq *
# constant cligr{1}
equation p2eq *
# constant cligr{1}
*
* Standard guess values aiming to make state 1 the low growth state.
*
compute v1=log(.7/.3)~~0.0
compute v2=log(.95/.05)~~0.0
compute mu(1)=-1.7,mu(2)=.3

boxjenk(ar=nlags) g
compute phi=||%beta(1),%beta(2),%beta(3),%beta(4)||
compute sigma=%seesq
*
maximize(parmset=tv+common,start=(pstar=TVPInit()),$
 method=bhhh,iters=20) logltvtp gstart gend

[TomDoan]

You need to do an @MSVARInitial at some point before you put in your own guess values.

[fan]

You need to do an @MSVARInitial at some point before you put in your own guess values.

Thank you for the advice. I have one more question would like to ask. Following the example of KIMNP078_1995.RPF, I am trying to modify the code to incorporate the possibility that the gap between average growth rats during expansion and recession could be smaller in more recent period. However, I am getting this error message

## SR10. Missing Values And/Or SMPL Options Leave No Usable Data Points.

I am not clear what caused this error as I do not see any missing values in the data set. Could you kindly take a look at my code? I am looking forward to hearing from you. Thank you in advance!!

Code: Select all

cal(m) 1948
open data filardo.dat
data(format=prn,org=cols) 1948:1 1992:8 year month ip cli dfi xli sp fed tspread
*
set ipgr = 100.0*log(ip/ip{1})
*
* The graph in the paper appears to be standardized (probably by dividing through
* by 2.0 times the sample standard deviation).
*
graph(footer="Figure 1. Log Growth Rate of Monthly IP")
# ipgr
*
stats(noprint) ipgr * 1959:12
compute stdearly=sqrt(%variance)
stats(noprint) ipgr 1960:1 *
compute stdlate=sqrt(%variance)
*
set ipgradjust = %if(t<=1959:12,ipgr*stdlate/stdearly,ipgr)
*
set g = ipgradjust
*
* Transform the information variables
*
set cligr   = 100.0*log(cli/cli{1})
set spgr    = 100.0*log(sp/sp{1})
set xlidiff = xli-xli{1}
set ffdiff  = fed-fed{1}
set tspdiff = tspread-tspread{1}
set dfismooth = dfi+2*dfi{1}+2*dfi{2}+dfi{3}
*
* And center them to zero means
*
dofor s = cligr spgr xlidiff ffdiff tspdiff dfismooth
   diff(center) s
end dofor s

********************************************************

* Common setup
*
compute nlags=4
*
source msvarsetup.src
compute gstart=1948:6,gend=1992:8
@MSVARSetup(lags=nlags)
# g
*********************************************************************
*Set up the shift dummy and the extra parameters on the shift
set dummy = t>=1985:4
dec vect mustar(2)
*
*************************************************************************
*Override the %MSVARFVec function. This returns the expanded set of likelihood
*functions across all state combinations.

function %MSVARFVec time
type vector %MSVARFVec
type integer time
*
local integer i j state
local vect u
*
dim %MSVARFVec(nexpand)
*
do i=1,nexpand
compute state=MSVARLagState(i,1)
compute u=%xt(y,time)-mu(state)(1)-mustar(state)*dummy(time)
do j=1,nlags
    compute statej=MSVARLagState(i,j+1), $
        u=u-phi(j)*(%xt(y,time-j)-mu(statej)(1)-mustar(statej)*dummy(time-j))
end do j
compute state=MSVARLagState(i,1)
compute %MSVARFVec(i)=exp(-.5*MSVARLogDetV(state)-.5*%qform(MSVARSigmaInv(state),u))
end do i
*
end

*********************************************************************
nonlin(parmset=common) mu phi sigma
@msvarinitial gstart gend

*********************************************************************
*
* Time varying probabilities
*
dec equation p1eq p2eq
dec vector   v1   v2
nonlin(parmset=tv) v1 v2
*********************************************************************
*
* This overrides the standard fixed definition
*
function  %MSVARPmat  time
type rect %MSVARPmat
type int  time
*
local integer i j
local rect    pexpand
local real    z
*
compute p(1,1)=%(z=exp(%dot(%eqnxvector(p1eq,time),v1)),z/(1+z))
compute p(1,2)=%(z=exp(%dot(%eqnxvector(p2eq,time),v2)),1/(1+z))
compute %MSVARInitTransition()
*
if nexpand==nstates {
   dim pexpand(nstates,nstates-1)
   ewise pexpand(i,j)=%if(i==nstates,-1.0,i==j)
   compute %MSVARPmat=pexpand*p
   ewise %MSVARPmat(i,j)=%if(i==nstates,1.0+%MSVARPmat(i,j),%MSVARPmat(i,j))
}
else {
   dim %MSVARPmat(nexpand,nexpand)
   ewise %MSVARPmat(i,j)=MSVARTransProbs(MSVARTransLookup(i,j))
}
end
*********************************************************************
*
* Initialize the probabilities using only the intercepts
*
function TVPInit
type vector TVPInit
*
local rect a
local integer i j
*
compute p=||%logistic(v1(1),1.0),1-%logistic(v2(1),1.0)||
compute TVPInit=%MSVARInit()
end
*********************************************************************
frml logltvtp = %MSVARPMat(t),log(%MSVarProb(t))
*********************************************************************
*
* Define the logistic index for the transitions
*
equation p1eq *
# constant cligr{1}
equation p2eq *
# constant cligr{1}
*
* Standard guess values aiming to make state 1 the low growth state.
*
compute v1=log(.7/.3)~~0.0
compute v2=log(.95/.05)~~0.0
compute mu(1)=-1.7,mu(2)=.3

maximize(parmset=tv+common,start=(pstar=TVPInit()),$
 method=bhhh,iters=20) logltvtp gstart gend

[TomDoan]

You never initialize MUSTAR and you don't include it in your parameter set.

[fan]

You never initialize MUSTAR and you don't include it in your parameter set.

Thank you so much for pointing out my silly mistakes. Have a nice day

[fan]

Dear Tom,

As you mentioned that to rescale that the growth rates of IP for period before 1960 would change the mean during that period, which wouldn't make much sense, could you please kindly suggest a better way to deal with data sets that have the similar pattern, dispersion in the early sample appears greater than in the later part? Many Thanks

[TomDoan]

There's no reason you can't have a variance which changes at a known location. It would require changing the @MSVARFVEC function.

[fan]

There's no reason you can't have a variance which changes at a known location. It would require changing the @MSVARFVEC function.

Thank you for your quick reply. Following your suggestion, I modified the @MSVARFVEC function. Could you please take a quick look and let me know whether it is correct? Many Thanks

Code: Select all

set dummy = t<=1959:12
dec real alfa
*********************************************************************
*
* Override the %MSVARFVec function. This returns the expanded set of
* likelihood functions across all state combinations.
*
function %MSVARFVec time
type vector  %MSVARFVec
type integer time
*
local integer i j state
local vect    u
*
dim %MSVARFVec(nexpand)
*
do i=1,nexpand
   compute state=MSVARLagState(i,1)
   compute u=%xt(y,time)-mu(state)(1)
   do j=1,nlags
      compute statej=MSVARLagState(i,j+1),$
        u=(1+alfa*time(t))*(u-phi(j)*(%xt(y,time-j)-mu(statej)(1))
   end do j
   compute state=MSVARLagState(i,1)
   compute %MSVARFVec(i)=exp(-.5*MSVARLogDetV(state)-.5*%qform(MSVARSigmaInv(state),u))
end do i
*
end

[TomDoan]

No. If you look at where you put the rescaling, it's applying the factor multiple times to the short lags. Plus, even if you get that corrected, it's not the same thing as the much simpler (and more justifiable) correction of the variances.

[fan]

No. If you look at where you put the rescaling, it's applying the factor multiple times to the short lags. Plus, even if you get that corrected, it's not the same thing as the much simpler (and more justifiable) correction of the variances.

Thank you for your quick response. I changed the code, but not sure whether the change is correct or not. Could you please kindly take a look? If the change is incorrect, could you please share a little a tip how to make the variance correction

Code: Select all

set dummy = t<=1959:12
dec real alfa
*********************************************************************
*
* Override the %MSVARFVec function. This returns the expanded set of
* likelihood functions across all state combinations.
*
function %MSVARFVec time
type vector  %MSVARFVec
type integer time
*
local integer i j state
local vect    u
*
dim %MSVARFVec(nexpand)
*
do i=1,nexpand
   compute state=MSVARLagState(i,1)
   compute u=%xt(y,time)-mu(state)(1)
   do j=1,nlags
      compute statej=MSVARLagState(i,j+1),$
        u=(u-phi(j)*(%xt(y,time-j)-mu(statej)(1))
   end do j
   compute state=MSVARLagState(i,1)
   compute %MSVARFVec(i)=exp(-.5*MSVARLogDetV(state)*(1+alfa*dummy(t))-.5*%qform(MSVARSigmaInv(state)(1+alfa*dummy(t))*,u))
end do i
*
end

[TomDoan]

What the correct change is depends upon whether 1+alfa*dummy(t) is supposed to inflate or deflate the variance. Think about what should happen to the log variance (first term) and the inverse variance (second term) when you multiply (or divide, depending upon the interpretation of the factor) the variance by a constant.

[fan]

What the correct change is depends upon whether 1+alfa*dummy(t) is supposed to inflate or deflate the variance. Think about what should happen to the log variance (first term) and the inverse variance (second term) when you multiply (or divide, depending upon the interpretation of the factor) the variance by a constant.

Thank you for sharing the tip. 1+alfa*dummy(t) is supposed to inflate the variance as the author points out that dispersion in the first third of the sample appears greater than in the later part. So, I need to inflation the log variance by multiplying 1+alfa*dummy(t), and deflating the inverse variance by dividing 1+alfa*dummy(t). Is that correct?

Code: Select all

 compute %MSVARFVec(i)=exp(-.5*MSVARLogDetV(state)*(1+alfa*dummy(t))-.5*%qform(MSVARSigmaInv(state)/(1+alfa*dummy(t)),u)

[TomDoan]

Multiplying the log by the inflation factor is *really* not the right thing to do. The first part of that is actually computing v^(-.5) (using exp(-.5*log(v)) so you either divide the whole factor by sqrt(1+alfa*dummy) or add log(1+alfa*dummy) to the v before it gets multiplied by the -.5.