**************************************************************
*This Program calculates a linear model with 1 latent variable
*for data with missing values (not all individuals are involved 
*in all trials).  The model is:
*            y_nt = a_t + b_t * X_n + e_nt
*  where:    y_nt is the observed variable for nth individual 
*                 tth trial (i.e. ADA score of member n in year t)
*            a_t =  the intercept (shift) for trial t
*            b_t =  the slope (stretch) for trial t
*            X_n =  the true unobserved (assumed constant) characteristic 
*                   of the nth individual in the sample.
*            e_nt = residual error assumes: e_nt~N(0,sigma) where sigma is
*                   is constant accross all individuals and trials.
* 
* In order to idenitify and estimate the model the following information
* is required:
*         o  An observed variable (e.g. ADA score in a given year)
*         o  A variable that uniquely identifies each individual  (e.g. MC)
*         o  A variable that uniquely identifies each trial (e.g. year chamber)
*         o  A scaling restriction setting a a_t=0 & b_t=1 for some trial t
*         
*  SYNTAX:
*     inflat <observed variable>, n(<n identifier>) t(<t identifier>)
*              i(<number of iterations>) cons(<t for a_t=0 b_t=1>) start 
*              disp
*
*     Options:  disp = display more info on what's happening in each interation
*               start = calculate start values for x 
*               i(<iterations>) = number of iterations, 100 is default
*
*
*
*  Jeff Lewis 
*  (jblewis@ucla.edu)
*  2000
* 
*************************************************************************


capture program drop inflat
program define inflat
version 5.0

set more 1

*************** Parse Arguements ***********************
local varlist "req"
local options "n(string) t(string) cons(string) disp start i(integer 100)" 
parse "`*'"
capture confirm `n'
capture confirm `t'
local year "`t'"
local member "`n'"

************** Setup temp vars & locals ****************
quietly {
capture drop a b
capture drop a_old b_old x_old
gen a_old = .
gen b_old = .
gen x_old = .
tempvar Y yr e2 sig mem isbase
local score "`varlist'"
if "`disp'" == "" { local disp "*" }
gen `sig' = .
gen `e2' = .
gen `Y' = .
}

************** Generate start values **********************
quietly {
if "`start'" != "" { egen x = mean(`score'), by(`member') }
gen a = 0
gen b = 1
}

*************  Give members a number from 1 to N **********
quietly {
sort `member' 
gen `mem' = 0
replace `mem' = 1 if `member' != `member'[_n-1]
replace `mem' = sum(`mem')
summ `mem'
local max = _result(6)
}


************  Give each year a number from 1 to T ******
sort `year'
quietly gen `yr'=0
quietly replace `yr' = 1 if `year'!=`year'[_n-1]
quietly replace `yr' = sum(`yr')
quietly gen `isbase' = `yr' if `year' == `cons'
quietly summ `isbase'
local isbasel = _result(6)
drop `isbase'
quietly summ `yr' 
local minyr = _result(5)
local maxyr = _result(6)

*********** Start the program running *******************
disp "Scaling: `score'"
disp "--------"
disp "Last t = `maxyr'"
disp "Last n = `max'"
disp "Base t = `isbasel'"
disp
pause

local itter = 1
while `itter' <= `i' {
  sort `year'
  disp
  disp "iteration = `itter'"
 ******* Get As & Bs ************
  local yrl = 1
  `disp' "    Finding As and Bs"
  `disp' "      Year: " _continue    
  while `yrl' <= `maxyr' {
    `disp' "`yrl'.." _continue
       if `yrl' != `isbasel' {
	 quietly {
	   reg `score' x if `yr'==`yrl'  
	   matrix b = get(_b)
	   replace a = b[1,2] if `yr' == `yrl' 
	   replace b = b[1,1] if `yr' == `yrl'
	 }
       }
  local yrl = `yrl' + 1
  }
  `disp'

  *******  Get Xs ******************************
  local meml = 1
  sort `mem' `yr'
  quietly replace `Y' = `score' - a
  `disp' "   Finding Xs"
  `disp' "     member: " _continue
  while `meml' < `max' {
     `disp' "`meml'.." _continue
     quietly {
	 capture reg `Y' b if `mem' == `meml', noc
	 if _rc == 0 {
         	matrix b = get(_b)
		replace x=b[1,1] if `mem' == `meml'
		}
	 else {
		`disp' in red " (1 obs) " _continue
		replace x = `Y'/b if `mem' == `meml'
		}
	 }
     local meml = `meml' + 1
   }
`disp'

*** generate llik *****
quietly replace `e2' = (`score' - a - b*x)^2
quietly replace `sig' = sum(`e2')
local sig2 = `sig'[_N] / _N
local llike = -_N * (ln(sqrt(2*_pi*sqrt(`sig2'))) + .5)

********* Setup for Next iteration *******
** Check for convergence **

local conv = .999999999
if `itter' > 1 {
	qui corr x x_old
	local corr = _result(4)
        disp in red "corr x = " _result(4)
	if `corr' > `conv' {
		qui corr b b_old
		local corr = _result(4)
		disp in red "corr b = " _result(4)
		if `corr' > `conv' {
			qui corr a a_old
			local corr = _result(4)
   		        disp in red "corr a = " _result(4)
                    	if `corr' > `conv' {
          			disp in green "Parameters converged"
   			        disp
			        local itter = `i' + 1
			        }
		        }
		}
	}			


** display end of iteration stuff **

disp
disp "Time = $S_TIME"
disp "Sigma = " %8.4f sqrt(`sig2')
disp "Log-likelihood = " %9.4f `llike'
disp in green "10th Obs:"   
disp in green "---------"
disp in green "Parameter" _column(10)"Value" _column(20)"Change"
disp in green "---------" _column(10)"-----" _column(20)"------"
disp in green "a" in yellow _column(10) a[10] _column(20) a[10] - a_old[10]
disp in green "b" in yellow _column(10) b[10] _column(20) b[10] - b_old[10]
disp in green "x" in yellow _column(10) x[10] _column(20) x[10] - x_old[10]
local itter = `itter' + 1
*quietly save c:\temp, replace
replace x_old = x
replace b_old = b
replace a_old = a
}

set more 0
drop x_old b_old a_old
end