////////////////////////////////////////////////////////////////////////////////
// STATA FOR  Hu, Y., Moffitt, R., & Sasaki, Y. (2019): Semiparametric
//            Estimation of the Canonical Permanent‐Transitory Model of Earnings
//            Dynamics. Quantitative Economics 10 (4), 1495-1536.
//
// Use it when you consider a state space model where the observed process is an
// sum of permanent and transitory components. The command draws the densities
// of the perment and transitory components .
////////////////////////////////////////////////////////////////////////////////
program define cdecompose, eclass
    version 14.2

    syntax varlist(min=2 numeric) [if] [in] [, p(real 1) q(real 1) delta(real 5) nboot(real 1000)]
    marksample touse

    gettoken depvar indepvars : varlist
    _fv_check_depvar `depvar'
    fvexpand `indepvars'
    local cnames `r(varlist)'

    tempname N b V rho mu kappa

	mata: estimate_moments("`depvar'", "`cnames'", ///
						   `p', `q', ///
						   `delta', `nboot', ///
				           "`touse'", "`N'", ///
						   "`b'", "`V'", "`rho'", "`mu'", "`kappa'")

	mata: bootstrap_moments("`depvar'", "`cnames'", ///
						    `p', `q', ///
						    `delta', `nboot', ///
				            "`touse'", "`N'", ///
						    "`b'", "`V'")

	local cnames U_Mean U_Std_Dev U_Skewness U_Kurtosis V_Mean V_Std_Dev V_Skewness V_Kurtosis
 	matrix colnames `b' = `cnames'
	matrix colnames `V' = `cnames'
	matrix rownames `V' = `cnames'



    ereturn post `b' `V', esample(`touse') buildfvinfo
    ereturn scalar N    = `N'
    ereturn scalar p    = `p'
    ereturn scalar q    = `q'
    ereturn local  cmd  "cdecompose"
	ereturn matrix rho  = `rho'

    ereturn display
end





mata:
////////////////////////////////////////////////////////////////////////////////
// Estimation
void estimate_moments( string scalar y1v,    		  string scalar y2v,
					   real scalar p,	  	  		  real scalar q,
					   real scalar delta,	  	  	  real scalar num_boot,
					   string scalar touse,   		  string scalar nname,
					   string scalar bname,   		  string scalar vname,
					   string scalar rhoname,         string scalar muname,
					   string scalar kappaname)
{
	printf("\n{hline 78}\n")
	printf("Executing: Hu, Y., Moffitt, R., & Sasaki, Y. (2019): Semiparametric Estimation\n")
	printf("           of the Canonical Permanent‐Transitory Model of Earnings Dynamics.  \n")
	printf("           Quantitative Economics 10 (4), 1495-1536.                          \n")
	printf("{hline 78}\n")

    y = st_data(., y1v, touse), st_data(., y2v, touse)
    n  = rows(y)
	// Set main index (t)
	tdx = p+q+1
	//tdx
	// Check the following is true
	if( cols(y) < p+q+2+q ){
		printf("Error: The number of variables must be at least p+2q+2=%f.\n", p+2*q+2)
	}
	// Demean y
	mean_y = mean(y[.,tdx])
	y = y :- mean_y
	// Set the 100 times the precision of numerical derivative in the frequency domain
	//delta = 100
	// Set slist
	slist = (-trunc(5*delta)..trunc(5*delta)) / 100
	slist_interval = slist[2] - slist[1]

	printf("Estimating the moments\n")
	////////////////////////////////////////////////////////////////////////////
	// Estimate rho
	if( p > 0 ){
		pprime = 0
		Delta0 = y[.,tdx+1-pprime]:-y[.,tdx-q]:-y[.,tdx-pprime]:+y[.,tdx-q-1]:+y[.,tdx-q]:-y[.,tdx-q-1]
		pprime = 1
		Delta = y[.,tdx+1-pprime]:-y[.,tdx-q]:-y[.,tdx-pprime]:+y[.,tdx-q-1]:+y[.,tdx-q]:-y[.,tdx-q-1]
		if( p > 1 ){
			for( pprime = 2 ; pprime <= p ; pprime++ ){
				 Delta = Delta, y[.,tdx+1-pprime]:-y[.,tdx-q]:-y[.,tdx-pprime]:+y[.,tdx-q-1]:+y[.,tdx-q]:-y[.,tdx-q-1]
			}
		}
		rho = luinv( (y[.,tdx :- ((q+1)..(p+q))])' * Delta :+ 0.05:*diag(J(p,1,n)) ) * (y[.,tdx :- ((q+1)..(p+q))])' * Delta0
	//	rho
	}

	////////////////////////////////////////////////////////////////////////////
	// Estimate mu
	mu = y[.,tdx+q+1] - y[.,tdx]
	if( p > 0 ){
		for( pprime = 1 ; pprime <= p ; pprime++ ){
			mu = mu :- rho[pprime] :* ( y[.,tdx+q+1-pprime] - y[.,tdx] )
		}
	}

	////////////////////////////////////////////////////////////////////////////
	// Estimate kappa
	kappa = -1
	if( p > 0 ){
		for( pprime = 1 ; pprime <= p ; pprime++ ){
			kappa = kappa + rho[pprime]
		}
	}
	//kappa

	////////////////////////////////////////////////////////////////////////////
	// Estimate phi_v
	phi_v_integrand = slist :* (0+0i)
	for( idx = 1 ; idx <= length(slist) ; idx++ ){
	s = slist[idx]
	phi_v_integrand[idx] = sum( (0+1i) :* mu :* exp( (0+1i):*s:*y[.,tdx] ) ) / sum( kappa :* exp( (0+1i):*s:*y[.,tdx] ) )
	}
	//phi_v_integrand'

	phi_v = slist :* (0+0i)
	middle = trunc(length(slist) / 2) + 1
	for( idx = middle + 1 ; idx <= length(slist) ; idx++ ){
		phi_v[idx] = phi_v[idx-1] + phi_v_integrand[idx] * slist_interval
	}
	for( idx = middle - 1 ; idx >= 1 ; idx-- ){
		phi_v[idx] = phi_v[idx+1] - phi_v_integrand[idx] * slist_interval
	}
	phi_v = exp( phi_v )
	//phi_v'

	////////////////////////////////////////////////////////////////////////////
	// Estimate phi_u
	phi_y = slist :* (0+0i)
	for( idx = 1 ; idx <= length(slist) ; idx++ ){
		s = slist[idx]
		phi_y[idx] = sum( exp( (0+1i):*s:*y[.,tdx] ) ) / n
	}
	phi_u = phi_y :/ phi_v
	//phi_u'

	////////////////////////////////////////////////////////////////////////////
	// Estimate derivatives of the characteristic functions
	phi_u_d1 = slist :* (0+0i)
	phi_v_d1 = slist :* (0+0i)
	phi_u_d2 = slist :* (0+0i)
	phi_v_d2 = slist :* (0+0i)
	phi_u_d3 = slist :* (0+0i)
	phi_v_d3 = slist :* (0+0i)
	phi_u_d4 = slist :* (0+0i)
	phi_v_d4 = slist :* (0+0i)

	for( idx = delta+1 ; idx <= length(slist)-delta ; idx++ ){
		phi_u_d1[idx] = (phi_u[idx+delta]-phi_u[idx-delta])/(slist[idx+delta]-slist[idx-delta])
		phi_v_d1[idx] = (phi_v[idx+delta]-phi_v[idx-delta])/(slist[idx+delta]-slist[idx-delta])
	}
	//phi_u_d1', phi_v_d1'

	for( idx = 2*delta+1 ; idx <= length(slist)-2*delta ; idx++ ){
		phi_u_d2[idx] = (phi_u_d1[idx+delta]-phi_u_d1[idx-delta])/(slist[idx+delta]-slist[idx-delta])
		phi_v_d2[idx] = (phi_v_d1[idx+delta]-phi_v_d1[idx-delta])/(slist[idx+delta]-slist[idx-delta])
	}
	//phi_u_d2', phi_v_d2'

	for( idx = 3*delta+1 ; idx <= length(slist)-3*delta ; idx++ ){
		phi_u_d3[idx] = (phi_u_d2[idx+delta]-phi_u_d2[idx-delta])/(slist[idx+delta]-slist[idx-delta])
		phi_v_d3[idx] = (phi_v_d2[idx+delta]-phi_v_d2[idx-delta])/(slist[idx+delta]-slist[idx-delta])
	}
	//phi_u_d3', phi_v_d3'

	for( idx = 4*delta+1 ; idx <= length(slist)-4*delta ; idx++ ){
		phi_u_d4[idx] = (phi_u_d3[idx+delta]-phi_u_d3[idx-delta])/(slist[idx+delta]-slist[idx-delta])
		phi_v_d4[idx] = (phi_v_d3[idx+delta]-phi_v_d3[idx-delta])/(slist[idx+delta]-slist[idx-delta])
	}
	//phi_u_d4', phi_v_d4'

	////////////////////////////////////////////////////////////////////////////
	// Estimate moments
	u_m1 = sum( phi_u_d1[((middle-delta)..(middle+delta))] ) / (0+1i) / (2*delta+1)
	u_m2 = sum( phi_u_d2[((middle-delta)..(middle+delta))] ) / (-1+0i) / (2*delta+1)
	u_m3 = sum( phi_u_d3[((middle-delta)..(middle+delta))] ) / (0-1i) / (2*delta+1)
	u_m4 = sum( phi_u_d4[((middle-delta)..(middle+delta))] ) / (1+0i) / (2*delta+1)

	v_m1 = sum( phi_v_d1[((middle-delta)..(middle+delta))] ) / (0+1i) / (2*delta+1)
	v_m2 = sum( phi_v_d2[((middle-delta)..(middle+delta))] ) / (-1+0i) / (2*delta+1)
	v_m3 = sum( phi_v_d3[((middle-delta)..(middle+delta))] ) / (0-1i) / (2*delta+1)
	v_m4 = sum( phi_v_d4[((middle-delta)..(middle+delta))] ) / (1+0i) / (2*delta+1)

	//u_m1, u_m2, u_m3, u_m4
	//v_m1, v_m2, v_m3, v_m4

	mean_u = u_m1
	sd_u = ( u_m2 - u_m1^2 )^0.5
	skew_u = ( u_m3 - 3*u_m2*mean_u + 3*u_m1*mean_u^2 - mean_u^3 ) / sd_u^3
	kurt_u = ( u_m4 - 4*u_m3*mean_u + 6*u_m2*mean_u^2 - 4*u_m1*mean_u^3 + mean_u^4 ) / sd_u^4

	mean_v = v_m1
	sd_v = ( v_m2 - v_m1^2 )^0.5
	skew_v = ( v_m3 - 3*v_m2*mean_v + 3*v_m1*mean_v^2 - mean_v^3 ) / sd_v^3
	kurt_v = ( v_m4 - 4*v_m3*mean_v + 6*v_m2*mean_v^2 - 4*v_m1*mean_v^3 + mean_v^4 ) / sd_v^4

	//mean_u, sd_u, skew_u, kurt_u
	//mean_v, sd_v, skew_v, kurt_v

	b = mean_u, abs(sd_u), skew_u, abs(kurt_u), mean_v, abs(sd_v), skew_v, abs(kurt_v)
	b = Re(b)

	////////////////////////////////////////////////////////////////////////////
	// Set
    st_numscalar(nname, n)
    st_matrix(bname, b)
	st_matrix(rhoname, rho)
	st_matrix(muname, mu)
    st_numscalar(kappaname, kappa)
}

////////////////////////////////////////////////////////////////////////////////
// Bootstrap
void bootstrap_moments( string scalar y1v,    		  string scalar y2v,
						real scalar p,	  	  		  real scalar q,
					    real scalar delta,	  	  	  real scalar num_boot,
						string scalar touse,   		  string scalar nname,
						string scalar bname,   		  string scalar vname)
{
    yy = st_data(., y1v, touse), st_data(., y2v, touse)
    n  = rows(yy)
	// Set main index (t)
	tdx = p+q+1
	//tdx
	// Demean y
	mean_yy = mean(yy[.,tdx])
	yy = yy :- mean_yy
	// Set the 100 times the precision of numerical derivative in the frequency domain
	//delta = 100
	// Set slist
	slist = (-trunc(5*delta)..trunc(5*delta)) / 100
	slist_interval = slist[2] - slist[1]

	printf("Estimating the variances of the moment estimates\n")
	//num_boot = 10
	blist = J(num_boot,8,0)
	for( boot_idx = 1 ; boot_idx <= num_boot ; boot_idx++ ){
		if( trunc(boot_idx/trunc(num_boot/10)) == boot_idx/trunc(num_boot/10) ){
			printf("                     %f%%\n",boot_idx/num_boot*100)
		}
		indices = trunc(uniform(n,1):*n):+1
		y = yy[indices,.]
		////////////////////////////////////////////////////////////////////////////
		// Estimate rho
		if( p > 0 ){
			pprime = 0
			Delta0 = y[.,tdx+1-pprime]:-y[.,tdx-q]:-y[.,tdx-pprime]:+y[.,tdx-q-1]:+y[.,tdx-q]:-y[.,tdx-q-1]
			pprime = 1
			Delta = y[.,tdx+1-pprime]:-y[.,tdx-q]:-y[.,tdx-pprime]:+y[.,tdx-q-1]:+y[.,tdx-q]:-y[.,tdx-q-1]
			if( p > 1 ){
				for( pprime = 2 ; pprime <= p ; pprime++ ){
					 Delta = Delta, y[.,tdx+1-pprime]:-y[.,tdx-q]:-y[.,tdx-pprime]:+y[.,tdx-q-1]:+y[.,tdx-q]:-y[.,tdx-q-1]
				}
			}
			rho = luinv( (y[.,tdx :- ((q+1)..(p+q))])' * Delta :+ 0.05:*diag(J(p,1,n)) ) * (y[.,tdx :- ((q+1)..(p+q))])' * Delta0
		//	rho
		}

		////////////////////////////////////////////////////////////////////////////
		// Estimate mu
		mu = y[.,tdx+q+1] - y[.,tdx]
		if( p > 0 ){
			for( pprime = 1 ; pprime <= p ; pprime++ ){
				mu = mu :- rho[pprime] :* ( y[.,tdx+q+1-pprime] - y[.,tdx] )
			}
		}

		////////////////////////////////////////////////////////////////////////////
		// Estimate kappa
		kappa = -1
		if( p > 0 ){
			for( pprime = 1 ; pprime <= p ; pprime++ ){
				kappa = kappa + rho[pprime]
			}
		}
		//kappa

		////////////////////////////////////////////////////////////////////////////
		// Estimate phi_v
		phi_v_integrand = slist :* (0+0i)
		for( idx = 1 ; idx <= length(slist) ; idx++ ){
		s = slist[idx]
		phi_v_integrand[idx] = sum( (0+1i) :* mu :* exp( (0+1i):*s:*y[.,tdx] ) ) / sum( kappa :* exp( (0+1i):*s:*y[.,tdx] ) )
		}
		//phi_v_integrand'

		phi_v = slist :* (0+0i)
		middle = trunc(length(slist) / 2) + 1
		for( idx = middle + 1 ; idx <= length(slist) ; idx++ ){
			phi_v[idx] = phi_v[idx-1] + phi_v_integrand[idx] * slist_interval
		}
		for( idx = middle - 1 ; idx >= 1 ; idx-- ){
			phi_v[idx] = phi_v[idx+1] - phi_v_integrand[idx] * slist_interval
		}
		phi_v = exp( phi_v )
		//phi_v'

		////////////////////////////////////////////////////////////////////////////
		// Estimate phi_u
		phi_y = slist :* (0+0i)
		for( idx = 1 ; idx <= length(slist) ; idx++ ){
			s = slist[idx]
			phi_y[idx] = sum( exp( (0+1i):*s:*y[.,tdx] ) ) / n
		}
		phi_u = phi_y :/ phi_v
		//phi_u'

		////////////////////////////////////////////////////////////////////////////
		// Estimate derivatives of the characteristic functions
		phi_u_d1 = slist :* (0+0i)
		phi_v_d1 = slist :* (0+0i)
		phi_u_d2 = slist :* (0+0i)
		phi_v_d2 = slist :* (0+0i)
		phi_u_d3 = slist :* (0+0i)
		phi_v_d3 = slist :* (0+0i)
		phi_u_d4 = slist :* (0+0i)
		phi_v_d4 = slist :* (0+0i)

		for( idx = delta+1 ; idx <= length(slist)-delta ; idx++ ){
			phi_u_d1[idx] = (phi_u[idx+delta]-phi_u[idx-delta])/(slist[idx+delta]-slist[idx-delta])
			phi_v_d1[idx] = (phi_v[idx+delta]-phi_v[idx-delta])/(slist[idx+delta]-slist[idx-delta])
		}
		//phi_u_d1', phi_v_d1'

		for( idx = 2*delta+1 ; idx <= length(slist)-2*delta ; idx++ ){
			phi_u_d2[idx] = (phi_u_d1[idx+delta]-phi_u_d1[idx-delta])/(slist[idx+delta]-slist[idx-delta])
			phi_v_d2[idx] = (phi_v_d1[idx+delta]-phi_v_d1[idx-delta])/(slist[idx+delta]-slist[idx-delta])
		}
		//phi_u_d2', phi_v_d2'

		for( idx = 3*delta+1 ; idx <= length(slist)-3*delta ; idx++ ){
			phi_u_d3[idx] = (phi_u_d2[idx+delta]-phi_u_d2[idx-delta])/(slist[idx+delta]-slist[idx-delta])
			phi_v_d3[idx] = (phi_v_d2[idx+delta]-phi_v_d2[idx-delta])/(slist[idx+delta]-slist[idx-delta])
		}
		//phi_u_d3', phi_v_d3'

		for( idx = 4*delta+1 ; idx <= length(slist)-4*delta ; idx++ ){
			phi_u_d4[idx] = (phi_u_d3[idx+delta]-phi_u_d3[idx-delta])/(slist[idx+delta]-slist[idx-delta])
			phi_v_d4[idx] = (phi_v_d3[idx+delta]-phi_v_d3[idx-delta])/(slist[idx+delta]-slist[idx-delta])
		}
		//phi_u_d4', phi_v_d4'

		////////////////////////////////////////////////////////////////////////////
		// Estimate moments
		u_m1 = sum( phi_u_d1[((middle-delta)..(middle+delta))] ) / (0+1i) / (2*delta+1)
		u_m2 = sum( phi_u_d2[((middle-delta)..(middle+delta))] ) / (-1+0i) / (2*delta+1)
		u_m3 = sum( phi_u_d3[((middle-delta)..(middle+delta))] ) / (0-1i) / (2*delta+1)
		u_m4 = sum( phi_u_d4[((middle-delta)..(middle+delta))] ) / (1+0i) / (2*delta+1)

		v_m1 = sum( phi_v_d1[((middle-delta)..(middle+delta))] ) / (0+1i) / (2*delta+1)
		v_m2 = sum( phi_v_d2[((middle-delta)..(middle+delta))] ) / (-1+0i) / (2*delta+1)
		v_m3 = sum( phi_v_d3[((middle-delta)..(middle+delta))] ) / (0-1i) / (2*delta+1)
		v_m4 = sum( phi_v_d4[((middle-delta)..(middle+delta))] ) / (1+0i) / (2*delta+1)

		//u_m1, u_m2, u_m3, u_m4
		//v_m1, v_m2, v_m3, v_m4

		mean_u = u_m1
		sd_u = ( u_m2 - u_m1^2 )^0.5
		skew_u = ( u_m3 - 3*u_m2*mean_u + 3*u_m1*mean_u^2 - mean_u^3 ) / sd_u^3
		kurt_u = ( u_m4 - 4*u_m3*mean_u + 6*u_m2*mean_u^2 - 4*u_m1*mean_u^3 + mean_u^4 ) / sd_u^4

		mean_v = v_m1
		sd_v = ( v_m2 - v_m1^2 )^0.5
		skew_v = ( v_m3 - 3*v_m2*mean_v + 3*v_m1*mean_v^2 - mean_v^3 ) / sd_v^3
		kurt_v = ( v_m4 - 4*v_m3*mean_v + 6*v_m2*mean_v^2 - 4*v_m1*mean_v^3 + mean_v^4 ) / sd_v^4

		//mean_u, sd_u, skew_u, kurt_u
		//mean_v, sd_v, skew_v, kurt_v

		b = mean_u, abs(sd_u), skew_u, abs(kurt_u), mean_v, abs(sd_v), skew_v, abs(kurt_v)
		b = Re(b)
		blist[boot_idx,.] = b
	}

	//sort(colshape(abs(blist),1),1)
	for( idx = 1 ; idx <= 8 ; idx++ ){
		large = blist[.,idx] :>= sort(blist[.,idx],1)[trunc(num_boot*0.995)+1]
		small = blist[.,idx] :<= sort(blist[.,idx],1)[trunc(num_boot*0.005)+1]
		blist[.,idx] = (!small :& !large) :* blist[.,idx] :+ small :* sort(blist[.,idx],1)[trunc(num_boot*0.005)+1] + large :* sort(blist[.,idx],1)[trunc(num_boot*0.995)+1]
	}
	V = ( ( blist' * blist ) :/ num_boot :- ( J(1,num_boot,1) * blist :/ num_boot)' * ( J(1,num_boot,1) * blist :/ num_boot) )

	////////////////////////////////////////////////////////////////////////////
	// Set
    st_numscalar(nname, n)
    st_matrix(vname, V)

	printf("\n{hline 78}\n")
	printf("ARMA(%f,%f) process - p=%f & q=%f                                         \n",p,q,p,q)
	printf("Number of cross-sectional observations:                               N=%f\n",n)
	printf("Number of time periods of input variables:                            T=%f\n",cols(yy))
	printf("Number of minimum time periods required for identification:           p+2q+2=%f\n",p+2*q+2)
	printf("{hline 78}\n")
	printf("Displayed are moments in the time period corresponding to the %f-th input var.\n",tdx)
	printf("U = Permanent Component\n")
	printf("V = Transitory Component\n")
}
end
////////////////////////////////////////////////////////////////////////////////