//////////////////////////////////////////////////////////////////////////////////////////////////////////////
/// Code version (v1.01): Feb 15, 2021 (v1.0: Feb 09,2021)												   ///
/// list of ado files: pcdid.ado, grtestsub.ado, pdd.ado												   ///
/// This code implements the PCDID estimator by Marc K. Chan and Simon Kwok,                               ///
/// "The PCDID Approach: Difference-in-Differences when Trends are Potentially Unparallel and Stochastic", ///
/// also previously circulated as "Policy Evaluation with Interactive Fixed Effects"                       ///
/// For more details, visit https://sites.google.com/site/marcchanecon/									   ///
//////////////////////////////////////////////////////////////////////////////////////////////////////////////

/***************************************************************************/
/* 20210207: prediction (postestimation: must run it after command pcdid) */
program define pdd
	version 14
	syntax namelist
	marksample touse
	disp as text "(WARNING: if you used an IF clause in estimation to select id's,     "
	disp as text "          the predictions generated by pdd will be unstable.)"
	if (e(cmd)=="pcdid") {
		cap drop `namelist'
		cap drop fproxy1-fproxy`e(factnum)'
		cap drop constant
		cap drop idnew
		cap drop treattmp
		gen constant = 1

		gen treattmp = `e(treatvar)'			/* generate temp treatment variable =0,0.5,1 (for appropriate sorting later) */
		if (e(treatlistnum) == 1) {				/* additional treated unit restriction by option */
			quietly replace treattmp = 0.5 if `e(treatvar)'==1 & (`e(treatlist)')
		}

		quietly merge n:1 `e(time)' using fproxy, keepusing (fproxy1-fproxy`e(factnum)')
		quietly drop _merge
		*sort `e(treatvar)' `e(id)' `e(time)'		/* important for alignment with matc and matb: sort by treatvar (0,1), then id, then time */
		*egen idnew = group(`e(treatvar)' `e(id)')
		sort treattmp `e(id)' `e(time)'				/* important for alignment with matc and matb: sort by treatvar (0,0.5,1), then id, then time */
		egen idnew = group(treattmp `e(id)')

		mkmat `e(time)' if idnew==1, matrix(timemat)	/* save time index for later merging */
		mat matall = e(matc) \ e(matb)				/* important: must be aligned with how the id's are arranged in the estimation command */

		forval i=1/`=scalar(e(Nc)+e(Ne))' {
			mkmat `e(indeps)' fproxy1-fproxy`e(factnum)' constant if idnew==`i', matrix(xmat)
			mat pdvec = xmat * matall[`i',1...]'
			if (`i'==1) {
				mat pdmat = pdvec
			}
			else{
				mat pdmat = pdmat, pdvec
			}
		}

		/* generate auxillary data set that contains predictions */
		preserve
		clear
		quietly svmat timemat, names(`e(time)')			/* store time index into data set */
		rename `e(time)'1 `e(time)'
		quietly svmat pdmat, names(`namelist')			/* store predicted values into data set */
		quietly reshape long `namelist', i(`e(time)') j(idnew)	/* reshape to long form for merging */
		quietly save yhatlong,replace
		restore

		/* merge main data with auxillary data set */
		quietly merge n:1 idnew `e(time)' using yhatlong
		drop _merge

		drop constant
		drop idnew
		drop treattmp			/* v1.01: refine over v1.0 */
		erase yhatlong.dta
	}
	else{
		disp "Error: must run pcdid command first"
	}
	disp as text "Predictions generated."
end program pdd