# st: RE: RE: Including a matrix in a monte carlo simulation

 From "Nevo, Dorit" To "'statalist@hsphsun2.harvard.edu'" Subject st: RE: RE: Including a matrix in a monte carlo simulation Date Wed, 16 Oct 2002 12:53:33 -0700

Thanks.

I'm not sure I'm going in the right direction so I'm pasting the full
program below. I am trying to calculate the mean and standard error of
sigma. I need quite a few calculations along the was so I tried to
generalize on the basic monte carlo example that I found in the reference
book and online.

By the time I create the matrix Fc should already be calculated. It works
fine without the matrix command so I'm not sure why the problem occurred.

Dorit

program define aes_monte
version 7.0
if "`1'" == "?" {
global S_1" mean var"
exit
}
drop _all
set obs 3000
gen B0 = -.0353521
gen R =  .5022559 + .0403609*invnorm(uniform())
gen Dc = -.3476753 + .0700479*invnorm(uniform())
gen Dk = .3394081 + .0281788*invnorm(uniform())
gen Bcc = -.0106232 + .0030351*invnorm(uniform())
gen Bkk = .0782757 +  .0023668*invnorm(uniform())
gen Bll = .0859557 + .0030879*invnorm(uniform())
gen Bck = .0217035 + .002353*invnorm(uniform())
gen Bcl = .0066442 +  .0040709*invnorm(uniform())
gen Bkl = -.1821175 + .0040324*invnorm(uniform())
gen C = 14.86249
gen lnC = 2.69884
gen K = 576.6429
gen lnK = 6.357223
gen L = 486.6226
gen lnL = 6.187489
gen V =exp(B0 - (1/R)*ln(Dc*(C^(-R)) + Dk*(K^(-R)) +
(1-Dc-Dk)*(L^(-R))) + Bcc*(lnC^2) + Bkk*(lnK^2) + Bll*(lnL^2) + Bck*lnC*lnK
+ Bcl*lnC*lnL + Bkl*lnK*lnL)

gen Z = Dc*C^(-R)+Dk*K^(-R)+(1-Dc-Dk)*L^(-R)
gen Fc = V*(Dc*C^(-R-1)/Z+Bck/C*lnK+Bcl/C*lnL+2*Bcc/C*lnC)
gen Fk = V*(Dk*K^(-R-1)/Z+Bck/K*lnC+Bkl/K*lnL+2*Bkk/K*lnK)
gen Fl = V*((1-Dc-Dk)*L^(-R-1)/Z+Bcl/L*lnC+Bkl/L*lnK+2*Bll/L*lnL)
gen Fcc = Fc^2/V - Fc/C +
V*((-R*Dc*C^(-R-2))/Z+(Dc^2*C^(-2*R-2)*R)/Z^2+2*Bcc/C^2)
gen Fkk = Fk^2/V - Fk/K +
V*((-R*Dk*K^(-R-2))/Z+(Dk^2*K^(-2*R-2)*R)/Z^2+2*Bkk/K^2)
gen Fll = Fl^2/V - Fl/L +
V*((-R*(1-Dc-Dk)*L^(-R-2))/Z+((1-Dc-Dk)^2*L^(-2*R-2)*R)/Z^2+2*Bll/L^2)
gen Fck = (Fc*Fk)/V + V*((R*Dc*Dk*C^(-R-1)*K^(-R-1)/Z^2)+Bck/(C*K))
gen Fcl = (Fc*Fl)/V +
V*((R*Dc*(1-Dc-Dk)*C^(-R-1)*L^(-R-1)/Z^2)+Bcl/(C*L))
gen Fkl = (Fk*Fl)/V +
V*((R*Dk*(1-Dc-Dk)*K^(-R-1)*L^(-R-1)/Z^2)+Bkl/(K*L))

matrix H = (0,Fc,Fk, Fl\Fc, Fcc,
Fck,Fcl\Fk,Fck,Fkk,Fkl\Fl,Fcl,Fkl,Fll)
gen detH = det(H)
matrix Hck = (0,Fc,Fl\Fk,Fkc,Fkl\Fl,Flc,Fll)
gen detHck = det(Hck)

gen sigma_ck =  ((C*Fc+K*Fk+L*Fl)/(C*K))*(-detHck/detH)
summarize sigma_ck
post `1' (r(mean)) (r(Var))
end

-----Original Message-----
From: Nick Cox [mailto:n.j.cox@durham.ac.uk]
Sent: Wednesday, October 16, 2002 12:40 PM
To: statalist@hsphsun2.harvard.edu
Subject: st: RE: Including a matrix in a monte carlo simulation

Nevo, Dorit
>
> I'm trying to calculate the determinant of a matrix within
> a Monte Carlo
> simulation program. The matrix takes variables that are
> calculated in
> previous steps in the program.
>
> The specific command I included is:
> 	matrix H =
> (0,Fc,Fk,Fl\Fc,Fcc,Fck,Fcl\Fk,Fck,Fkk,Fkl\Fl,Fcl,Fkl,Fll)
> 	gen detH = det(H)
>
> Fc, Fk, etc are calculated earlier in the program file.
>

Stata won't allow you to insert a whole variable
in this position, even if it only obtains one
observation. The strange message is emitted because
Stata is prepared to put named scalars in like
this, but evidently you have no scalar named -Fc-.

You can specify Fc[1] or whatever to insert
particular values for a variable.

On the whole it seems likely that it would be much
better style and more efficient for the upstream
calculation to produce scalars, not variables, but
in the absence of most of your code, that's
a guess.

But you _are_ putting a determinant -- one number --
in the new variable -detH-, and that is not
necessary on the face of it.

What might make sense to loop around roughly like
this

gen detH = .

forval i = 1 / whatever {

do some stuff

form the matrix H

qui replace detH = det(H) in `i'

}

Nick
n.j.cox@durham.ac.uk
*
*   For searches and help try:
*   http://www.stata.com/support/faqs/res/findit.html
*   http://www.stata.com/support/statalist/faq
*   http://www.ats.ucla.edu/stat/stata/
*
*   For searches and help try:
*   http://www.stata.com/support/faqs/res/findit.html
*   http://www.stata.com/support/statalist/faq
*   http://www.ats.ucla.edu/stat/stata/