I think I've managed to run the code completely. However, I have a problem with the last point where the impulse response is configured, I can not generate the multiple graph, which is the last code.
I suspect that my problem is with "LOWER", but I can not identify what this is.
I give you the code and sample data I'm working with. This works properly and runs in its entirety.
I also leave you, again, my code and my data. Please help me to know the error, I think I'm very close to being able to complete this exercise.
Code: Select all
OPEN DATA "Z:\Users\Icaro\Dropbox\Trabajo de Grado MaestriÌa\replication_IER\indexsa.xls"
CALENDAR 1980 1 4
ALL 2000:04
DATA(FORMAT=XLS,ORG=COLUMNS) 1980:01 2000:04 USINF JPINF GEINF UKINF FRINF ITINF CAINF CAINFR GEINFR FRINFR $
ITINFR JPINFR UKINFR USINFR CAINFW GEINFW FRINFW ITINFW JPINFW UKINFW USINFW CAEG GEEG FREG ITEG JPEG $
UKEG USEG CAMG GEMG FRMG ITMG JPMG UKMG USMG CAREE GEREE FRREE ITREE JPREE UKREE USREE CAGDPG GEGDPG $
FRGDPG ITGDPG JPGDPG UKGDPG USGDPG CASLOPE GESLOPE FRSLOPE ITSLOPE JPSLOPE UKSLOPE USSLOPE CAGAP GEGAP $
FRGAP ITGAP JPGAP USGAP UKGAP UKSTG CPG USINFR1 esree esinf esinfr eseg eseg1 esgdpg nlree nlinf nlinfr $
nlgdpg nleg EA_gdpg EA_hicpg EA_pcdg EA_eg EA_stn
********
* Se hacen los ajustes de la muestra de los datos
seed 11
cal 1980 1 4
all 2002:4
*Inicialmente este comando ALL estaba de la forma all 14 2002:4, lo cual me desconfiguraba la selección de la muestra. Por tanto se modificó eliminando el número 14 que seguía el comando ALL
source(noecho) cndition.src
compute beg=1980:1
compute tempend=2000:4
compute dat=tempend-beg
compute obs=tempend-beg
compute nstep=16
smpl beg tempend
*******
*Se definen los parámetros a utlizar en el modelo. Inicia con el parámetro alter1, el cual hace referencia a una dummy asociada con la varianza. si alter1=1 se estima el sigma^2 desde los datos. Si alter1=0 se fija el sigma^2=0.0001
*Se define dim_lam igual a 2 y se cambia a 1/1 si lamda tiene la mitad 1/2 de componentes comunes. Es decir, exógenas que afectan a las unidades (países)
*Se define zeta=.00 para definir el escenario no informativo en sigma(t)
*Se define exact como una dummy, que será 1 si la factorización es exacta y será igual a 0 si no lo es.
compute alter1=1
compute dim_lam=2
compute zeta=.00
compute exact=1
******
*Ahora lo que hacemos establecer las condiciones geneales del modelo,
*países=7, variables=4 y dimensiones=variables*países. También se definen los min, med y max de los
*rezagos, tanto en las variables endógenas como en las exógenas.
compute n=7
compute g=4
compute ng=n*g
compute minlag=0
compute medlag=1
compute maxlag=1
compute lag_exo=1
compute nlags=p=maxlag-medlag+1
compute p1=lag_exo-medlag+1
compute nlags1=3
compute d1=0
compute d=d1*p1+1
compute k1=n*g*p+d
compute ndraws=2100
compute ndraws1=30
compute mindraws=2000
compute begdraws=ndraws-mindraws
compute k_1=n+g+dim_lam
display k_1
******
*Se definen los valores de algunos parámetros
compute scala1=0.000001
compute scala2=0.000001
compute scala3=0.000001
compute arch3=1.0
compute arch4=0.0
compute w_0=n*g+50
compute v_0=1000000, d_0=1.0
******
*Datos
open data indexsa.xls
data(org=obs,for=xls)
close data
******
*Reorganizamos algunas series para reducir afectaciones al ejercicio
*La tasa de crecimiento de Alemania presenta un pico muy fuerte
*entre los años 1990 y 1993, por tanto se procede a "fraccionarla"
statistics gegdpg beg tempend
statistics gegdpg beg 1990:4
compute a1=%mean
statistics gegdpg 1992:1 2000:4
compute a2=%mean
statistics geeg beg 1990:4
compute a3=%mean
statistics geeg 1992:1 2000:4
compute a4=%mean
set gegdpg 1991:1 1991:2 =a1
set gegdpg 1991:3 1991:4 =a2
set geeg 1991:1 1991:2 =a3
set geeg 1991:3 1991:4 =a4
statistics gegdpg beg tempend
*Se fija una inflación promedio de los siete países igual a cero?
*Todo el siguiente segmento está en comentario, así que por ahora no se tendrá
*en consideración
/*
set inf_averg7=0.0
do i=1,n
set inf_averg7=inf_averg7+(1/n)*([series] i)
end
*Estandarización de las variables
dofor i = usinf usinfr useg usgdpg jpinf jpinfr jpeg jpgdpg $
geinf geinfr geeg gegdpg ukinf ukinfr ukeg ukgdpg $
frinf frinfr freg frgdpg itinf itinfr iteg itgdpg $
cainf cainfr caeg cagdpg cpg inf_averg7
stat(noprint) ([series] i)
set ([series] i) = (([series] i)-%mean)/sqrt(%variance)
*set ([series] i) = ([series] i)/sqrt(%variance)
end dofor
*/
******
*Estimamos ahora sigma^2 y Q_0
* Q_0 es sigma_eps y omega también se incia como Q_0
**** Ahora se estima un var que va desde 1 hasta n*g=7*4=28.
*se usan las 28 variables, un rezago para cada una de ellas y se incluye además la constante.
*Posterioemente se realiza una estimación para toda la muestra, y de esta se extrae la matriz
*de varianzas y covarianzas sigma_eps. Luego decalramos una matriz omega, de simensión (ng,ng)
*y se designa como la inversa de sigma_eps, es decir, la inversa de la matriz se var y covar
*Luego sacamos una matriz (29,28) con los coeficientes obtenidos para todas las regresoras
*del VAR
system(model=var1) 1 to n*g
variables usinf usinfr useg usgdpg jpinf jpinfr jpeg jpgdpg $
geinf geinfr geeg gegdpg ukinf ukinfr ukeg ukgdpg $
frinf frinfr freg frgdpg itinf itinfr iteg itgdpg $
cainf cainfr caeg cagdpg
lags 1 to 1
det constant
end(system)
estimate(noprint,noftest,outsigma=sigma_eps) beg tempend
display sigma_eps
dec symm omega(ng,ng)
compute omega = inv(sigma_eps)
display omega
compute beta_ols1=%modelgetcoeffs(var1)
display beta_ols1
display %rows(beta_ols1) %cols(beta_ols1)
*******
*Ahora sigma^2
*No sé de dónde es alter. Se supone que yo lo fijo, cómo lo reconoce
*el programa a partir de lo que ya hemos calculado.
*Lo anterior ya no es importante, puesto que se plantea el modelo de forma condicional
*para tener en consideración las dos opciones del valor de alter= 0 o 1. Si es 1, sigma se
*calcula a partir de los datos, de acuerdo con la fórmula que allí se muesta. Si es 0, entonces
*se e asigna el valor sigma=0.0001. %seesq= error estandar de la estimación al cuadrado.
compute s_e=0.0
if alter1==1
{
dofor i = usinf usinfr useg usgdpg jpinf jpinfr jpeg jpgdpg $
geinf geinfr geeg gegdpg ukinf ukinfr ukeg ukgdpg $
frinf frinfr freg frgdpg itinf itinfr iteg itgdpg $
cainf cainfr caeg cagdpg
linreg(noprint) i
# constant i{1 to nlags1}
compute s_e = s_e+(1./(n*g))*(sqrt(%seesq))
end dofor
}
if alter1==0
compute s_e = 0.0001
end if
*
compute k_1=n+g+dim_lam
dec symm big_phi(k_1,k_1) big_phiinv(k_1,k_1)
comp big_phi = %const(.0)
*En los siguientes pasos, se adecuan las dimensiones de las matrices, sus composiciones y
*sus respectivas inversas. Esto con el propósito de poder llenar toda la matriz Phi mayúscula.
*Los phi_1, phi_2 y phi_3 formarán la diagonal de la matriz big_phi
dec symm phi_1(dim_lam,dim_lam)
comp phi_1 = %mscalar(scala1)
compu phi_1inv = inv(phi_1)
dec symm phi_2(g,g) phi_2inv(g,g)
dec symm phi_3(n,n) phi_3inv(n,n)
compute phi_2 = %mscalar(scala2)
comp phi_2inv = inv(phi_2)
compute phi_3 = %mscalar(scala3)
comp phi_3inv = inv(phi_3)
do i = 1,dim_lam
comp big_phi(i,i) = phi_1(i,i)
end
display big_phi
do i = 1,g
comp big_phi(i+dim_lam,i+dim_lam) = phi_2(i,i)
end
display big_phi
do i = 1,n
comp big_phi(i+g+dim_lam,i+g+dim_lam) = phi_3(i,i)
end
display big_phi
comp big_phiinv = inv(big_phi)
**************
*Se configuran las matrices para la estimación preliminar de OLS para inicializar los índices
*y para usar en en MCMC
**************
smpl(reglist)
# usinf{minlag to maxlag} usinfr{minlag to maxlag} useg{minlag to maxlag} usgdpg{minlag to maxlag} $
jpinf{minlag to maxlag} jpinfr{minlag to maxlag} jpeg{minlag to maxlag} jpgdpg{minlag to maxlag}$
geinf{minlag to maxlag} geinfr{minlag to maxlag} geeg{minlag to maxlag} gegdpg{minlag to maxlag} $
ukinf{minlag to maxlag} ukinfr{minlag to maxlag} ukeg{minlag to maxlag} ukgdpg{minlag to maxlag} $
frinf{minlag to maxlag} frinfr{minlag to maxlag} freg{minlag to maxlag} frgdpg{minlag to maxlag} $
itinf{minlag to maxlag} itinfr{minlag to maxlag} iteg{minlag to maxlag} itgdpg{minlag to maxlag} $
cainf{minlag to maxlag} cainfr{minlag to maxlag} caeg{minlag to maxlag} cagdpg{minlag to maxlag} $
constant
MAKE(transpose) Y
# usinf usinfr useg usgdpg $
jpinf jpinfr jpeg jpgdpg geinf geinfr geeg gegdpg $
ukinf ukinfr ukeg ukgdpg frinf frinfr freg frgdpg $
itinf itinfr iteg itgdpg cainf cainfr caeg cagdpg
comp ti = %cols(y)
*El siguiente comando en mayúscula indica la creación de un vector, cuyas entradas,
*a su vez, son matrices rectagulares. Para los siguientes, son las definiciones de los
*vectores y sus dimensiones.
DEC VEC[REC] XX(N+1)
dec vec hetero_old(ti)
dec vec hetero_can(ti)
dec vec hetero(ti)
dec rec hetero_dra(ti,ndraws)
comp hetero = %const(s_e)
*display hetero
dec vec s_t(ti)
dec vec[vec] ss_t(i)
*
dec rec[series] general1(ti,dim_lam)
dec rect[series] general2(ti,n)
dec rec[series] variables1(ti,g)
dec rec[series] vari_count1(ti,ng)
*
*Tanto la X mayúscula como la x minúscula son requeridas para correr el modelo. La X small
*está compuesta de vectores y esta a su vez es de dimensión 83x1. La X grande(big) está
*compuesta de rectangulares, y también es de dimensión 83x1
DEC VECT[vec] SMALLX(ti)
display %rows(smallx) %cols(smallx)
DEC VECT[RECT] BIGX(TI)
***Lo que sigue a continuación se requiere para no generar incompatibilidades de dimensiones
*not interdep/interdep
*Construimos la matriz X con los rezagos de las variables endógenas.
make X
# usinf{medlag to maxlag} usinfr{medlag to maxlag} useg{medlag to maxlag} usgdpg{medlag to maxlag} $
jpinf{medlag to maxlag} jpinfr{medlag to maxlag} jpeg{medlag to maxlag} jpgdpg{medlag to maxlag}$
geinf{medlag to maxlag} geinfr{medlag to maxlag} geeg{medlag to maxlag} gegdpg{medlag to maxlag} $
ukinf{medlag to maxlag} ukinfr{medlag to maxlag} ukeg{medlag to maxlag} ukgdpg{medlag to maxlag} $
frinf{medlag to maxlag} frinfr{medlag to maxlag} freg{medlag to maxlag} frgdpg{medlag to maxlag} $
itinf{medlag to maxlag} itinfr{medlag to maxlag} iteg{medlag to maxlag} itgdpg{medlag to maxlag} $
cainf{medlag to maxlag} cainfr{medlag to maxlag} caeg{medlag to maxlag} cagdpg{medlag to maxlag} $
constant
*display x
**Ahora calculamos el número de observaciones que se pueden usar y el número de variables
*por cada ecuación
compu ti = %rows(x)
comp k = %cols(x)
display %rows(x) %cols(x)
*Lo que sigue indica que para cada uno de las observaciones (83), la dimensión de smallx será
*será de 29 (las 28 variables más la constante). Los siguientes dos DO se indican para poder
*aplicar el modelo a cualquier ejercicio, indepentientemente de cuáles sean los n,g y k
DO T = 1,TI
DIM SMALLX(T)(K)
END DO T
DO T = 1,TI
DO J = 1,K
COMPUTE SMALLX(T)(J)=x(T,J)
END DO J
END DO T
display smallx
*Ahora se configura la matriz X(t-s), es decir la matriz de los rezagos,
*lo cual se necesita para la reparametrización del modelo
*El número de rezagos dependerá del valor que se le de a s
do t = 1,ti
compute s_t(t) = 1.+s_e*%dot(smallx(t),smallx(t))
display s_t(t)
end
DO T = 1,T
DIM BIGX(T)(N*G,N*G*K)
COMPUTE BIGX(T) = %KRONEKER(%IDENTITY(N*G),TR(SMALLX(T)))
*display %rows(bigx(t)) %cols(bigx(t)) bigx(t)
END DO T
*Ahora se configura el vector de las variables endógenas, de acuerdo
*con las dimesiones de ti y k
DEC VECT[VECT] SMALLY(TI)
do t= 1,ti
overlay y(1,t) with SMALLY(t)(N*G)
end do t
**********************************************
*LAMNDA, RHO Y ALPHA*
**********************************************
*Para la configuarión del VAR como un sistema de índices, debemos organizar lamnda, roh y alpha
*Para eso definimos sus dimensiones y componentes, para después poder multiplicarlos por X
*Es importante tener presente que estas tres matrices nos permitirán como establecer
*Los efectos que queremos analizar, pues serán las matrices de ceros y unos que nos permitirán
*tener encuenta o no algunas variables.
*****MATRIZ LAMNDA
dec rec xsi_11(k,dim_lam) xsi_1(n*g*k,dim_lam)
compute xsi_1 = %const(.0), xsi_11 = %const(.0)
*
dec vec e_p(nlags)
do ll = 1,nlags
compute e_p(ll) = ll**(-0.5)
display e_p(ll) %rows(e_p) %cols(e_p)
end
dec vec e_p1(p1)
do ll = 1,p1
compute e_p1(ll) = ll**(-0.5)
display e_p1(ll) %rows(e_p1) %cols(e_p1)
end
dec vec e_d(d1) e_ng(n*g) e_ngp(n*g*p+d1*p1) xsi11(k) xsi1_1(n*g*k)
compute e_ng = %const(1.), e_d = %const(1.)
display e_ng %rows(e_ng) %cols(e_ng) e_d %rows(e_d) %cols(e_d)
if dim_lam==1
{
compute xsi11(k) = 1.
comp engp = %kroneker(e_ng,e_p)
display engp %rows(engp) %cols(engp)
comp edp = %kroneker(e_d,e_p1)
display edp %rows(edp) %cols(edp)
do j = 1,n*g*p
comp e_ngp(j) = engp(j,1)
display e_ngp(j)
end
do j = 1,d1*p1
comp e_ngp(n*g*p+j) = edp(j,1)
display e_ngp(n*g*p+j) %rows(e_ngp) %cols(e_ngp)
end
do j = 1,n*g*p+d1*p1
compute xsi11(j) = e_ngp(j)
display xsi11(j) %rows(xsi11) %cols(xsi11)
end
comp xsi1_1 = %kroneker(e_ng,xsi11)
display xsi1_1 %rows(xsi1_1) %cols(xsi1_1)
do j = 1,N*g*k
compute xsi_1(j,1) = xsi1_1(j)
display xsi_1(j,1) %rows(xsi_1) %cols(xsi_1)
end
}
*
***Ahora se configura escenario de Europa VS el resto del mundo
*Es decir, cuando la dimensión de lamnda (componentes comunes-predeterminados) es igual a 2
*
if dim_lam==2
{
do j = 1,d1*p1
comp xsi_11(n*g*p+j,1) = xsi_11(n*g*p+j,2) = e_p1(j)
end
comp xsi_11(k,1) = xsi_11(k,2) = 1.
do j = 1,g
do ll = 1,p
comp xsi_11((j-1)*p+ll,1) = e_p(ll)
end
end
do j = 1,g
do ll = 1,p
comp xsi_11(g*p+(j-1)*p+ll,1) = e_p(ll)
end
end
do j = 1,g
do ll = 1,p
comp xsi_11(2*g*p+(j-1)*p+ll,2) = e_p(ll)
end
end
do j = 1,g
do ll = 1,p
comp xsi_11(3*g*p+(j-1)*p+ll,1) = e_p(ll)
end
end
do j = 1,g
do ll = 1,p
comp xsi_11(4*g*p+(j-1)*p+ll,2) = e_p(ll)
end
end
do j = 1,g
do ll = 1,p
comp xsi_11(5*g*p+(j-1)*p+ll,2) = e_p(ll)
end
end
do j = 1,g
do ll = 1,p
comp xsi_11(6*g*p+(j-1)*p+ll,1) = e_p(ll)
end
end
do j = 1,g
do jj = 1,k
comp xsi_1((j-1)*k+jj,1) = xsi_11(jj,1)
end
end
do j = 1,g
do jj = 1,k
comp xsi_1(g*k+(j-1)*k+jj,1) = xsi_11(jj,1)
end
end
do j = 1,g
do jj = 1,k
comp xsi_1(3*g*k+(j-1)*k+jj,1) = xsi_11(jj,1)
end
end
do j = 1,g
do jj = 1,k
comp xsi_1(6*g*k+(j-1)*k+jj,1) = xsi_11(jj,1)
end
end
do j = 1,g
do jj = 1,k
comp xsi_1(2*g*k+(j-1)*k+jj,2) = xsi_11(jj,2)
end
end
do j = 1,g
do jj = 1,k
comp xsi_1(4*g*k+(j-1)*k+jj,2) = xsi_11(jj,2)
end
end
do j = 1,g
do jj = 1,k
comp xsi_1(5*g*k+(j-1)*k+jj,2) = xsi_11(jj,2)
end
end
}
end if
*display xsi_1 %rows(xsi_1) %cols(xsi_1)
****MATRIZ RHO
dec rec xsi2(g*k,g) xsi_2(n*g*k,g)
comp xsi2 = %const(.0), xsi_2 = %const(.0)
*display %rows(xsi2) %cols(xsi2) xsi2 %rows(xsi_2) %cols(xsi_2) xsi_2
dec vec e_p(nlags)
do ll = 1,nlags
compute e_p(ll) = ll**(-0.5)
end
*display e_p
dec vec[vec] a_xsi2(g)
do i = 1,g
dim a_xsi2(i)(k)
comp a_xsi2(i) = %const(.0)
do j = 1,d1*p1
comp a_xsi2(i)(n*g*p+j) = e_p1(j)
end
comp a_xsi2(i)(k) = 1.
end
*display %rows(a_xsi2(1)) %cols(a_xsi2) a_xsi2
do j = 1,g
do i = 1,n
do ll = 1,p
comp a_xsi2(j)((j-1)*p+(i-1)*g*p+ll) = e_p(ll)
end
end
end
*display %rows(a_xsi2) %cols(a_xsi2) a_xsi
do j = 1,g
do i = 1,k
comp xsi2((j-1)*k+i,j) = a_xsi2(j)(i)
end
end
*display %rows(xsi2) %cols(xsi2) xsi2
dec vec e_n(n)
compute e_n = %const(1.)
*display e_n
comp xsi_2 = %kroneker(e_n,xsi2)
*display %rows(xsi_2) %cols(xsi_2)
****MATRIZ ALPHA
dec rec xsi_3(n*g*k,n)
comp xsi_3 = %const(.0)
display %rows(xsi_3) %cols(xsi_3)
dec vec e_p(nlags)
do ll = 1,nlags
compute e_p(ll) = ll**(-0.5)
end
dec vec e_g(g) e_pg(p*g)
compute e_g = %const(1.)
*display e_g
comp e_pg = %kroneker(e_g,e_p)
*dislplay %rows(e_pg) %cols(e_pg) e_pg
dec vec[vec] a_xsi3(n)
do i = 1,n
dim a_xsi3(i)(k)
comp a_xsi3(i) = %const(.0)
do j = 1,d1*p1
comp a_xsi3(i)(n*g*p+j) = e_p1(j)
end
comp a_xsi3(i)(k) = 1.
end
*dislplay %rows(a_xsi3) %cols(a_xsi3)
do i = 1,n
do j = 1,g*p
comp a_xsi3(i)((i-1)*g*p+j) = e_pg(j)
end
end
dislplay %rows(a_xsi3) %cols(a_xsi3) a_xsi3
dec vec[vec] a_a_xsi3(n)
do i = 1,n
dim a_a_xsi3(i)(g*k)
comp a_a_xsi3(i) = %kroneker(e_g,a_xsi3(i))
end
do i = 1,n
do j = 1,g*k
comp xsi_3((i-1)*g*k+j,i) = a_a_xsi3(i)(j)
end
end
*display %rows(xsi_3) %cols(xsi_3) xsi_3
**********************************************
*MATRIZ DE REGRESORES
**********************************************
*El siguiente comando nos permite crear la matriz de regresores que se construirá a partir de
*un vector, cuyas entradas serán vectores rectangulares. Creamos tres de estos vectores,
*cada uno de dimensión ti=83, que son el número de observaciones por cada variable
*Luego se calculan los diferentes productos entre la matriz de regresores y los lamnda, rho
* y alpha, los cuales presentarán los efectos que se quieren analizar.
DEC VECT[RECT] BIG_W(TI) BIG_Z(TI) BIG_A(TI)
DO T = 1,ti
DIM BIG_W(T)(N*G,dim_lam) BIG_Z(T)(N*G,G) BIG_A(T)(N*G,N)
comp big_w(t) = bigx(t)*xsi_1
comp big_z(t) = bigx(t)*xsi_2
comp big_a(t) = bigx(t)*xsi_3
end
*display %rows(bigx(t)) %cols(bigx(t)) %rows(xsi_1) %cols(xsi_1)
****MATRIZ M_1
*Es decir, la reconfiguración del modelo en su versión más reducida
*Recordemos que esta matriz contiene los "coeficientes" que determinan los índices, en realidad
*en un vector compuesto por xsi_1, xsi_2 y xsi_3, para explicar a b(t)
dec rect m_1(n*g*k,n+g+dim_lam)
comp m_1 = %const(.0)
*display %rows(m_1) %cols(m_1) m_1
do i = 1,n*g*k
do lam = 1,dim_lam
compute m_1(i,lam) =xsi_1(i,lam)
end do lam
end do i
do i = 1,ng*k
do j = dim_lam+1,g+dim_lam
compute m_1(i,j) = xsi_2(i,j-dim_lam)
end do j
end do i
do i = 1,ng*k
do ii = g+dim_lam+1,dim_lam+g+n
compute m_1(i,ii) = xsi_3(i,ii-g-dim_lam)
end do ii
end do i
*
dec vect[rect] X_M(ti)
DO T = 1,ti
DIM x_m(T)(N*G,dim_lam+g+n)
comp X_M(t) = %const(.0)
do i = 1,n*g
do lam = 1,dim_lam
compute X_M(t)(i,lam) =big_w(t)(i,lam)
end
end do i
do i = 1,ng
do j = dim_lam+1,g+dim_lam
compute X_M(t)(i,j) = big_z(t)(i,j-dim_lam)
end do j
end do i
do i = 1,ng
do ii = g+dim_lam+1,dim_lam+g+n
compute X_M(t)(i,ii) = big_a(t)(i,ii-g-dim_lam)
end do i
end do i
end do t
display %rows(x_m) %cols(x_m)
display %rows(x_m(83))
*display x_m(83)
**********************************************************
*AQUÍ EMPIEZA LA CONFIGURACIÓN PARA MONTECARLO
**********************************************************
*Estimación de los índices usando ols secuenciales. (en prueba de muestras)
*Usa la estimación para inicializar MCMC
*Para calcular esto, basta tomar el modelo en su expresión extensa y despejar cada letra griega
*desde una forma matricial para obtener esta forma funcional.
dec vect[vect] lambda(ti) rho(ti) alfa_1(ti)
do t = 1,ti
dim lambda(t)(dim_lam) rho(t)(g) alfa_1(t)(n)
compute lambda(t) = inv(tr(big_w(t))*big_w(t))*tr(big_w(t))*smally(t)
compute eps1 = smally(t)-big_w(t)*lambda(t)
compute rho(t) = inv(tr(big_z(t))*big_z(t))*tr(big_z(t))*eps1
compute eps2 = eps1-big_z(t)*rho(t)
compute alfa_1(t) = inv(tr(big_a(t))*big_a(t))*tr(big_a(t))*eps2
end do t
comp t1 = ti
*
dec vect[vect] theta_ka(t1) theta_mc(t1)
do t = 1,t1
DIM theta_ka(t)(k_1)
dim theta_mc(t)(k_1)
COMPUTE THETA_KA(T) = %const(.0), THETA_MC(T) = %CONST(.0)
do la = 1,dim_lam
compute theta_ka(t)(la) = lambda(t)(la)
compute theta_mc(t)(la) = lambda(t)(la)
end
do j = 1,g
COMPUTE theta_ka(t)(j+dim_lam) = rho(t)(j)
compute theta_mc(t)(j+dim_lam) = rho(t)(j)
end do j
do i = 1,n
COMPUTE theta_ka(t)(g+dim_lam+i) = alfa_1(t)(i)
compute theta_mc(t)(g+dim_lam+i) = alfa_1(t)(i)
end do j
end do t
*display %rows(theta_mc)
*******
*Ahora se configuran las matrices necesarias para hacer el muestreo de Gibbs
******
dec vect[vect] u2(ti)
dec rect part5(n*g,n*g)
dec real part6 part7 part8 part9
compute def1 = t1+w_0
dec symm ss0hatinv(n*g,n*g)
*
* Se configuran las matrices para los valores futuros de los indicadores
*(opuesto a rezagos)
*
dec rec[vec] lamdraw(ti,ndraws)
dec rec[vec] rhodraw(ti,ndraws)
dec rec[vec] alfadraw(ti,ndraws)
dec vec phi1dra(ndraws) phi2dra(ndraws) phi3dra(ndraws)
dec vect[symm] omega_ka(t1) omega_mc(t1)
dec vect[rect] emme_mc(t1)
do t = 1,t1
dim omega_ka(t)(k_1,k_1) omega_mc(t)(k_1,k_1)
dim emme_mc(t)(k_1,k_1)
compute omega_ka(t)=%const(0.0), omega_mc(t)=%const(0.0)
compute emme_mc(t)=%const(0.0)
end do t
dec vec loglike_dr(ndraws)
dec vec[symm] omega_dra(ndraws)
*********************************
*COMIENZO RUTINA DE MCMC
*********************************
*
comp hetero = %const(1.)
do draws = 1,ndraws
disp draws
do t = 1,ti
dim lamdraw(t,draws)(dim_lam)
dim rhodraw(t,draws)(g)
dim alfadraw(t,draws)(n)
dim u2(t)(k_1)
end do t
* do draws1 = 1,ndraws1
*disp draws draws1
*
** kalman filter
*
compute omega_ka(1) = %IDENTITY(K_1)
compute omegaTT_1 = omega_ka(1)+big_phi
compute emme_mc(1) = omega_ka(1)*inv(omegaTT_1)
compute omega_mc(1) = omega_ka(1)-%mqform(omegaTT_1,tr(emme_mc(1)))
comp u2(1) = %ran(1.)
comp loglike_kf = 0.0
do t = 2,t1
compute delta_t = arch4**t+arch3*(1./(1.-arch4))*(1.-arch4**t)
compute omegaTT_1 = omega_ka(t-1)+delta_t*big_phi
compute effeTT_1 = X_M(t)*omegaTT_1*tr(X_M(t))+(1./hetero(t))*inv(omega)
compute kappaTT_1 = omegaTT_1*tr(X_M(t))*inv(effeTT_1)
compute imenoka = %identity(k_1)-kappaTT_1*X_M(t)
compute omega_ka(t) = imenoka*omegaTT_1
compute resi = smally(t)-X_M(t)*theta_ka(t-1)
compute theta_ka(t) = theta_ka(t-1)+kappaTT_1*resi
compute delta_t1 = arch4**(t+1)+arch3*(1./(1.-arch4))*(1.-arch4**(t+1))
compute omegaTT_11 = omega_ka(t)+delta_t1*big_phi
compute emme_mc(t) = omega_ka(t)*inv(omegaTT_11)
compute omega_mc(t)= omega_ka(t)-%mqform(omegaTT_11,tr(emme_mc(t)))
compute resi1 = smally(t)-X_M(t)*theta_ka(t)
comp loglike_temp = -0.5*n*g*log(2.*%pi)-0.5*log(%det(effeTT_1))-0.5*%qform(inv(effeTT_1),resi1)
comp loglike_kf = loglike_kf+loglike_temp
release resi resi1 omegaTT_1 effeTT_1 kappaTT_1 imenoka omegatt_11
comp u2(t) = %ran(1.)
end do t
comp loglike_dr(draws) = loglike_kf
do draws1 = 1,ndraws1
*DISP DRAWS DRAWS1
comp u2(t1) = %ran(1.)
compute theta_mc(t1) = theta_ka(t1) ;*+%decomp(omega_ka(t1))*u2(t1)
do sp = 1,DIM_LAM
compute lambda(t1)(sp) = theta_mc(t1)(sp)
end do sp
do sp = dim_lam+1,g+dim_lam
compute rho(t1)(sp-dim_lam) = theta_mc(t1)(sp)
end do sp
do sp = g+dim_lam+1,k_1
compute alfa_1(t1)(sp-g-dim_lam) = theta_mc(t1)(sp)
end do sp
do t = 2,t1
comp u2(t1-t+1) = %ran(1.)
compute themean=theta_ka(t1-t+1)+emme_mc(t1-t+1)* $
(theta_mc(t1-t+2)-theta_ka(t1-t+1))
compute theta_mc(t1-t+1) = themean+%decomp(omega_mc(t1-t+1))*u2(t1-t+1)
compute themean = %const(.0)
*compute theta_mc(t1-t+1) = theta_ka(t1-t+1) ;*+%decomp(omega_ka(t1))*u2(t1)
do sp = 1,dim_lam
compute lambda(t1-t+1)(sp) = theta_mc(t1-t+1)(sp)
end do sp
do sp = dim_lam+1,g+dim_lam
compute rho(t1-t+1)(sp-dim_lam) = theta_mc(t1-t+1)(sp)
end do sp
do sp = g+dim_lam+1,k_1
compute alfa_1(t1-t+1)(sp-g-dim_lam) = theta_mc(t1-t+1)(sp)
end do sp
end do t
* compute theta_ka(1) = theta_mc(1)
compute part5 = %const(.0)
compute part6 = .0
compute part7 = .0, part8 = .0, part9 = 0.0
comp sum_sig = 0.0, sum_eome = 0.0
DO T = 2,t1
compute eps1 = smally(t)-big_w(t)*lambda(t)-big_z(t)*rho(t)-big_a(t)*alfa_1(t)
comp eps_eps = eps1*tr(eps1)
comp part6 = zeta*s_t(t)+%qform(omega,eps1)
* compute eps1 = %const(.0)
compute p00 = 2./(part6)
if exact ==0
{
compute hetero(t) = p00*%rangamma(.5*(ng+zeta))
}
if exact ==1
{
compute hetero(t) = 1.0
}
compute part5 = part5+(hetero(t))*eps_eps
compute sum_sig = sum_sig+log(1./hetero(t))
comp sum_eome = sum_eome +(hetero(t))*%qform(omega,eps1)
compute delta_t = arch4**t+arch3*(1./(1.-arch4))*(1.-arch4**t)
compute resth = (1./delta_t)*(theta_mc(t)-theta_mc(t-1))
compute resth1 = (1./delta_t)*(lambda(t)-lambda(t-1))
compute resth2 = (1./delta_t)*(rho(t)-rho(t-1))
compute resth3 = (1./delta_t)*(alfa_1(t)-alfa_1(t-1))
compute part7 = part7+%dot(resth1,resth1)
compute part8 = part8+%dot(resth2,resth2)
compute part9 = part9+%dot(resth3,resth3)
end do t
compute seps1 = part7+d_0
compute p11= 2./(seps1)
compute phi1 = p11*%rangamma(.5*(t1+v_0))
compute phi_1 = %mscalar(1./phi1)
compute seps2 = part8+d_0
compute p2= 2./(seps2)
compute phi2 = p2*%rangamma(.5*(t1*g+v_0))
compute phi_2 = %mscalar(1./phi2)
compute seps3 = part9+d_0
compute p3= 2./(seps3)
compute phi3 = p3*%rangamma(.5*(t1*n+v_0))
compute phi_3 = %mscalar(1./phi3)
comp big_phi = %const(.0)
do i = 1,dim_lam
comp big_phi(i,i) = phi_1(i,i)
end do lam
do i = 1,g
comp big_phi(i+dim_lam,i+dim_lam) = phi_2(i,i)
end
do i = 1,n
comp big_phi(i+g+dim_lam,i+g+dim_lam) = phi_3(i,i)
end
comp big_phiinv = inv(big_phi)
compute werre_0 = w_0*sigma_eps
compute ss0hatinv = werre_0+part5
compute Swish = tr(%decomp(inv(ss0hatinv)))
compute omega = %mqform(%RANWISHART(N*g,def1),Swish)
end do draws1
comp omega_dra(draws) = omega
compute phi1dra(draws) = 1./phi1
compute phi2dra(draws) = 1./phi2
compute phi3dra(draws) = 1./phi3
do t = 2,ti
compute lamdraw(t,draws) = lambda(t)
compute rhodraw(t,draws) = rho(t)
compute alfadraw(t,draws) = alfa_1(t)
end do t
end do draws
******************************************
****** FIN DE LA RUTINA MCMC ***************
*****************************************
*Incertudumbre procedente de MCMC
comp nvar = ng
smpl beg tempend
dec vec[vec] theta_all(ndraws) delta_all(ndraws)
do draws = 1,ndraws
dim theta_all(draws)(k_1)
dim delta_all(draws)(n*g*k)
do sp = 1,dim_lam
compute theta_all(draws)(sp) = lamdraw(t1,draws)(sp)
end
do sp = dim_lam+1,dim_lam+g
compute theta_all(draws)(sp) = rhodraw(t1,draws)(sp-dim_lam)
end
do sp = g+dim_lam+1,k_1
compute theta_all(draws)(sp) = alfadraw(t1,draws)(sp-g-dim_lam)
end
comp delta_all(draws) = m_1*theta_all(draws)
end
*
system(model=varmodel1) 1 to nvar
variables usinf usinfr useg usgdpg jpinf jpinfr jpeg jpgdpg $
geinf geinfr geeg gegdpg ukinf ukinfr ukeg ukgdpg $
frinf frinfr freg frgdpg itinf itinfr iteg itgdpg $
cainf cainfr caeg cagdpg
lags 1 to nlags
det constant
end(system)
declare vect[rect] responses(ndraws)
declare rect[series] impulses(nvar,nvar)
dec rec betadraw(k,ng)
infobox(action=define,progress,lower=1,upper=ndraws) 'impulse responses'
do draws = 1,ndraws
do j = 1,ng
* comp delta_all(draws)((j-1)*(k+nlags)+1) = 0.6
do i = 1,k
compute betadraw(i,j) = delta_all(draws)((j-1)*k+i)
end
end
compute %modelsetcoeffs(varmodel1,(betadraw))
*
** Cálculo del bloque de Choleski
comp omega_11 = %xsubmat(inv(omega_dra(draws)),1,g,1,g)
comp omega_22 = %xsubmat(inv(omega_dra(draws)),1+g,n*g,1+g,n*g)
comp omega_12 = %xsubmat(inv(omega_dra(draws)),1,g,1+g,n*g)
comp omega_21 = %xsubmat(inv(omega_dra(draws)),1+g,n*g,1,g)
comp E_E = omega_22-omega_21*inv(omega_11)*omega_12
comp eigdec_om = %EIGDECOMP(omega_11)
comp eigdec_e = %EIGDECOMP(E_E)
comp lambda_om = %diag(eigdec_om(1))
comp lambda_e = %diag(eigdec_e(1))
comp sqrt_om = eigdec_om(2)*%sqrt(lambda_om)*tr(eigdec_om(2))
comp sqrt_e = eigdec_e(2)*%sqrt(lambda_e)*tr(eigdec_e(2))
comp fact = omega_21*inv(omega_11)*sqrt_om
comp chole = %blockglue(|| %BLOCKGLUE(|| sqrt_om , %zeros(g,(n-1)*g) ||) | %blockglue(|| fact, sqrt_e ||) ||)
** Finaliza el cálculo del bloque de Choleski
comp sigm = inv(omega_dra(draws))
display sigm
*display responses
source condition.src
*También se puede calcular choleski de esta forma
@cndition(decomp=chole) varmodel1 nstep tempend+1 sigm 2 condfore
# usgdpg tempend+1 usgdpg(tempend)+1.
* # usgdpg tempend+2 usgdpg(tempend)+.5
* # usgdpg tempend+3 usgdpg(tempend)+.25
* # usgdpg tempend+4 usgdpg(tempend)
# useg tempend+1 useg(tempend)+1.
* # useg tempend+2 usgdpg(tempend)+.5
forecast(model=varmodel1,results=forecasts) nvar nstep tempend+1
/*
set prova tempend+1 tempend+nstep = condfore(24)-forecasts(24)
set prova1 tempend+1 tempend+nstep = condfore(4)-forecasts(4)
graph(key=below) 2 ;#condfore(24) tempend+1 tempend+nstep ;#forecasts(24) tempend+1 tempend+nstep
graph(key=below) 2 ;#prova1 tempend+1 tempend+nstep ;#prova tempend+1 tempend+nstep
*/
* Almacenar los impulso respuesta como la diferencia entre el pronóstico condicional
*e incondicional
dim responses(draws)(nvar,nstep)
do i = 1,nvar
do j = 1,nstep
comp responses(draws)(i,j)=condfore(i)(tempend+j)-forecasts(i)(tempend+j)
end
end
infobox(current=draws)
end do draws
infobox(action=remove)
** compute here IRF and bands for each variable
** 4 = US gdp; 8 = JP gdp; 12 = DE gdp ; 16 = UK; 20 = FR gdp; 24 = IT gdp; 28 = CA gdp
dec vect[series] upper(nvar) lower(nvar) resp(nvar)
comp j = 16
do j = 1,nvar
smpl 1 ndraws
compute minlower=maxupper=0.0
clear lower(j) upper(j) resp(j)
do i=1,nstep
set work 1 ndraws = responses(t)(j,i)
compute frac=%fractiles(work,||.05,0.5,.95||)
compute lower(j)(i)=frac(1)
compute upper(j)(i)=frac(3)
compute resp(j)(i)=frac(2)
end do i
compute maxupper=%max(maxupper,%maxvalue(upper(j)))
compute minlower=%min(minlower,%minvalue(lower(j)))
/*
smpl 1 nstep
graph(ticks,min=minlower,max=maxupper,number=0) 3 j i
# resp(j)
# upper(j) / 2
# lower(j) / 2
*/
end
/*
do j = 4,28,4
print / resp(j) upper(j) lower(j)
end
*smpl 1 nstep
graph 2
# resp(4)
# resp(8)
# resp(12)
# resp(16)
# resp(20)
# resp(24)
# resp(28)
*/
COMP MINLOWER = -.5, MAXUPPER = 1.
smpl 1 nstep
spgraph(hfields=3,vfields=3,header='Responses of GDP growth to a shock to real US variables ')
graph(ticks,min=minlower,max=maxupper,number=0,header='US') 3 2 1
# resp(4)
# upper(4) / 2
# lower(4) / 2
graph(ticks,min=minlower,max=maxupper,number=0,header='Japan') 3 1 2
# resp(8)
# upper(8) / 2
# lower(8) / 2
graph(ticks,min=minlower,max=maxupper,number=0,header='Germany') 3 1 3
# resp(12)
# upper(12) / 2
# lower(12) / 2
graph(ticks,min=minlower,max=maxupper,number=0,header='UK') 3 2 2
# resp(16)
# upper(16) / 2
# lower(16) / 2
graph(ticks,min=minlower,max=maxupper,number=0,header='France') 3 2 3
# resp(20)
# upper(20) / 2
# lower(20) / 2
graph(ticks,min=minlower,max=maxupper,number=0,header='Canada') 3 3 2
# resp(28)
# upper(28) / 2
# lower(28) / 2
graph(ticks,min=minlower,max=maxupper,number=0,header='Italy') 3 3 3
# resp(24)
# upper(24) / 2
# lower(24) / 2
spgraph(done)