*! version 6.0.0  15Apr2000        STB-58: dm82
capture program drop simuped
program define simuped2
	version 6
	qui {

		gettoken age1 0 : 0, parse(" ,")
        	gettoken std1 0 : 0, parse(" ,")
        	gettoken age2 0 : 0, parse(" ,")
        	gettoken std2 0 : 0, parse(" ,")

        	if (`age1'<0 | `age2'<0 | `std1'<0 | `std2'<0) {
                	di in red "negative numbers invalid"
                	exit 498
        	}

		syntax [, Saving(string) Reps(int 100) Alle(real 0.1) Sib(real 3)]  

		if "`saving'"~="" {
			parse "`saving'", parse (" ")
			local  output= "`1'"
		}	
		else {
			local  output= "temp.dta"
		}
		
		tempname nsimu p q simu

		if (`reps' < 1) {
			di in red "reps() required"
			exit 198
		}
		scalar `nsimu' = `reps'

		if (`sib' < 0) {
			di in red "Sibship size negative"
			exit 198
		}

		if (`alle' < 0) {
			di in red "Allele frequency negative"
			exit 198
		}
		scalar `p' = `alle'
		scalar `q' = 1-`p'

		local nob=int(10*`sib'+2)
		set obs `nob'
		gen mu = `sib'
		
		tempvar simu degree id female nochild x g y age  i

		postfile simuped2 famid id degree female g age using `output', replace

		gen `simu'=0
		while (`simu'<`nsimu') {		/* begin simu */

		 	replace `simu'=`simu'+1
	
			tempvar degree id female nochild x g y age  i
			
	
			/*********** SIZE OF SIBSHIP ************/
			rndpoix mu
			gen byte `nochild'=xp 
			replace `nochild'=. if _n~=1 		/* denoted by 1st person */
			drop xp


			gen byte `degree'=1 if _n<=2
			replace `degree'=2 if _n>2 & _n<=2+`nochild'[1]
			gen byte `id'=_n if _n<=2+`nochild'[1]


					
			/************* SEX AND AGE *****************/
			gen byte `female'=0 if `degree'==1
			replace `female'=1 if _n==2 & `degree'==1		/* father first */
			
			gen `x'=uniform() if `degree'==2	
			replace `female'=0 if `x'<0.5 & `degree'==2
			replace `female'=1 if `x'>=0.5 & `degree'==2

			gen int `age'=max(0,int(`age1'+`std1'*invnorm(uniform())+0.5)) if `degree'==1
			replace `age'=max(0,int(`age2'+`std2'*invnorm(uniform())+0.5)) if `degree'==2

			/************* GENERATION 1: GENOTYPE ************/
			replace `x'=uniform() if `degree'==1
			gen byte `g'=11 if `x'<=`p'*`p' & `degree'==1
			replace `g'=12 if `x'>`p'*`p' & `x'<=`p'*`p'+2*`p'*`q' & `degree'==1
			replace `g'=22 if `x'>`p'*`p'+2*`p'*`q' & `degree'==1
				

			/************* GENERATION 2: GENOTYPE ************/


			replace `x'=uniform() if `degree'==2
			gen `y'=uniform() if `degree'==2

			
			#delimit ;
			replace `g'=11 if `degree'==2 &
						(`g'[1]==11 | `g'[1]==12 & `x'<0.5) &
                        		        (`g'[2]==11 | `g'[2]==12 & `y'<0.5);
			replace `g'=22 if `degree'==2 &
						(`g'[1]==22 | `g'[1]==12 & `x'>=0.5) & 
                          		        (`g'[2]==22 | `g'[2]==12 & `y'>=0.5);
			replace `g'=12 if `degree'==2 &
						(`g'[1]==11 | `g'[1]==12 & `x'<0.5) & 
                        		        (`g'[2]==22 | `g'[2]==12 & `y'>=0.5);
			replace `g'=12 if `degree'==2 &
						(`g'[1]==22 | `g'[1]==12 & `x'>=0.5) &
                          	  	        (`g'[2]==11 | `g'[2]==12 & `y'<0.5);
						
			#delimit cr


			/************ OUTPUT TO STATA FILE ***********/
			gen `i'=1			/* the ith observation */

			while (`i'<=2+`nochild'[1]) {
	
				post simuped2 `simu'[`i'] `id'[`i'] `degree'[`i'] /*
					*/ `female'[`i'] `g'[`i'] `age'[`i']
				replace `i'=`i'+1
			}
			noi di in green `simu' _skip(1) _con		/* display dots during simu */


		drop `degree' `id' `female' `nochild' `x' `g' `y' `age' `i'


		} 		/* end SIMU */

		postclose simuped2

		use `output', clear
		gen str2 genotype="AA" if g==11
		replace genotype="Aa" if g==12
		replace genotype="aa" if g==22
		drop g
		save `output', replace
		drop _all
	}
	end