精华内容
下载资源
问答
  • matlab对于probit模型怎么求统计值(赠probit程序)
    2021-04-20 01:53:51

    matlab对于probit模型怎么求统计值,如pseudo-R2。。。(附带程序来自网络)

    function [estimator,cov_Hessian,ME1,ME2,ME_std] = PROBIT(Y,X,method_flag)

    % Purpose:

    % Estimate Probit model and the marginal effects

    % -----------------------------------

    % Model:

    % Yi* = Xi * Beta + ui , where normalized ui ~ N(0,1)

    % Yi* is unobservable.

    % If Yi* > 0, we observe Yi = 1; If Yi* <= 0, we observe Yi = 0

    % -----------------------------------

    % Algorithm:

    % Maximum likelihood with analytic gradients

    % -----------------------------------

    % Usage:

    % Y = dependent variable (n * 1 vector)

    % X = regressors (n * k matrix)

    % method_flag = numeric gradients (1), analytic gradients (2, default), analytic Hessian (3)

    % -----------------------------------

    % Returns:

    % estimator = estimator corresponding to the k regressors

    % cov_Hessian = covariance matrix of the estimator

    % ME1 = marginal effects (average data)

    % ME2 = marginal effects (individual average)

    % ME_std = standard error of ME1

    % -----------------------------------

    % Notes:

    % Probit model is subject to normalization.

    % The variance of disturbances is set to 1, and a constant is added to X.

    %

    % Written by Hang Qian, Iowa State University

    % Contact me:  matlabist@gmail.com

    if nargin < 2;  error('Incomplete data');   end;

    if nargin < 3;  method_flag = 2;    end

    try

    JustTest = normcdf(1);

    catch

    disp(' ')

    disp('Oooopse, Matlab statistics toolbox is not installed.')

    disp('You may download a compatibility package on my website.')

    disp('http://www.public.iastate.edu/~hqi/toolkit')

    error('Program exits.')

    end

    [nrow_x,ncol_x] = size(X);

    [nrow_y,ncol_y] = size(Y);

    if nrow_x < ncol_x;    X = X';    ncol_x = nrow_x;end

    if nrow_y < ncol_y;    Y = Y';    ncol_y = nrow_y;end

    if ncol_x < ncol_y;    Y_temp = Y;    Y = X;    X = Y_temp;end

    [nobs,nreg] = size(X);

    Y_unique = unique(Y);

    if length(Y_unique) ~= 2

    disp('Ooopse, dependent variable should be binary.');

    disp('For illustration purpose, Y are grouped to binary data by its mean value')

    end

    Y = (Y > mean(Y));

    add_const = 0;

    const_check = all(X==1);

    if ~any(const_check)

    disp('A constant terms is added to X due to normalization. ')

    X = [ones(nobs,1),X];

    nreg = nreg + 1;

    const_check = [true,const_check];

    add_const = 1;

    end

    c_initial = (X'*X)\(X'*Y)*(rand-0.5)*10;

    switch method_flag

    case 1

    options = optimset('LargeScale','off','MaxFunEvals',1000,'Display','off');

    case 2

    options = optimset('LargeScale','off','GradObj','on','MaxFunEvals',1000,'Display','off','DerivativeCheck','off');

    case 3

    options = optimset('LargeScale','on','GradObj','on','Hessian','on','MaxFunEvals',2000,'MaxIter',1000,'Display','off','DerivativeCheck','off');

    end

    try

    [estimator,log_like] = fminunc(@(c)ML_PROBIT(c,Y,X,method_flag),c_initial,options);

    catch

    method_flag = 1;

    [estimator,log_like] = fminsearch(@(c)ML_PROBIT(c,Y,X,method_flag),c_initial);

    end

    q = 2*Y - 1;

    lambda = q .* normpdf(q.*(X*estimator)) ./ normcdf(q.*(X*estimator));

    Hessian = -transpose(lambda.*(lambda+X*estimator)*ones(1,nreg).*X)*X;

    fisher_info = -Hessian;

    cov_Hessian = inv(fisher_info);

    % Gradient_indiv = lambda*ones(1,nreg).*X;

    % cov_BHHH = inv(Gradient_indiv'*Gradient_indiv);

    std_c = sqrt(diag(cov_Hessian));

    t_stat = estimator./std_c;

    X_bar = mean(X);

    ME1 = normpdf(X_bar*estimator) * estimator;

    ME2 = mean(normpdf(X*estimator)) * estimator;

    Jocobian = normpdf(X_bar*estimator) * (eye(nreg)-X_bar*estimator*(estimator*X_bar));

    ME_var = Jocobian * cov_Hessian * Jocobian';

    ME_std = sqrt(diag(ME_var));

    ME1(const_check) = NaN;

    ME2(const_check) = NaN;

    ME_std(const_check) = NaN;

    eval([char([81 72 49 61]),'[87 114 105 116 116 101 110 32 98 121];'])

    eval([char([81 72 50 61]),'[32 72 97 110 103 32 81 105 97 110];'])

    result = cell(nreg+1,7);

    result(1,:) = {' ','Estimator','SE','z-stat','ME(avg. data)','ME_std','ME(ind. avg.)'};

    for m=1:nreg

    result(m+1,1) = {['C(',num2str(m-add_const),')']};

    result(m+1,2:7) = {estimator(m),std_c(m),t_stat(m),ME1(m),ME_std(m),ME2(m)};

    end

    disp(' ')

    disp(result)

    dummy_check = all((X==1) | (X==0)) & (~const_check);

    for m = 1:nreg

    if dummy_check(m)

    disp(['C(',num2str(m-add_const),') corresponds to binary data.'])

    X_dummy1 = X;

    X_dummy0 = X;

    X_dummy1(:,m) = ones(nobs,1);

    X_dummy0(:,m) = zeros(nobs,1);

    Y_margin_modify_data_bar = normcdf(mean(X_dummy1)*estimator) - normcdf(mean(X_dummy0)*estimator);

    Y_margin_modify_indi_bar = mean(normcdf(X_dummy1*estimator) - normcdf(X_dummy0*estimator));

    disp(['Adjusted marginal effects (average data) is ',num2str(Y_margin_modify_data_bar)])

    disp(['Adjusted marginal effects (individual average) is ',num2str(Y_margin_modify_indi_bar)])

    disp(' ')

    end

    end

    disp(['Log likelihood: ',num2str(-log_like)])

    disp(' ')

    fwrite(1, char([QH1,QH2,10,13]))

    end

    %-----------------------------------

    %  ML_PROBIT

    %-----------------------------------

    function [log_like,Gradient_c,Hessian_c] = ML_PROBIT(c,Y,X,method_flag)

    q = 2*Y-1;

    probit_F = normcdf(q.*(X*c));

    log_like = sum(log(probit_F));

    log_like = -log_like;

    if method_flag >= 2

    lambda = q.*normpdf(q.*(X*c))./probit_F;

    Gradient_c = lambda'*X;

    Gradient_c = -Gradient_c;

    end

    if method_flag == 3

    [nobs,nreg] = size(X);

    Hessian_c = transpose(lambda.*(lambda+X*c)*ones(1,nreg).*X)*X;

    Hessian_c = -Hessian_c;

    end

    end

    更多相关内容
  • matlab做probit回归代码RIS 时空概率 Philipp Hunziker 2018 年 8 月 12 日 这个包实现了 (FHC) 中引入的时空自回归概率模型。 估计基于递归重要性采样 (RIS) 最大模拟似然 (MSL) 程序。 该代码松散地基于 Franzese ...
  • 1、广义线性回归 广义线性模型有三个组成部分: (1) 随机部分, 即变量所属的指数族分布 族成员, 诸如正态分布, 二项分布, ...经过Probit变换和Logit变换,两种模型可以写成: 多变量情况: logit回归 probit回归 3
  • 二元选择(Probit&Logit)模型.doc
  • 为了评估二元经济结果之间的因果关系,我们考虑对面板数据的双变量动态概率模型进行估计,该模型具有特殊性,可以说明动态过程的初始条件。 由于似然函数的难解形式是二维积分,因此我们使用一种近似方法:自适应...
  • 针对矿井底板突水预测精确度要求高的特点,提出了二值Probit回归矿井突水预测模型。依据概率论的相关知识构建模型参数的极大似然估计法,并通过遗传算法求模型参数的最优解。利用所建模型对建模样本的突水性进行反向...
  • 建立了交通事故受伤人数预测的Ordered Probit模型,应用极大似然方法进行了模型标定.应用模型分析了受伤人数的影响因素,计算了各影响因素的边际贡献,并进行了受伤人数的预测.结果表明,对受伤人数影响较大的因素有路表...
  • 主要讨论如何采用 Probit模型作为研究方法进行实证研究。对 Probit模型的估计方法进行了详细解释,并对其估计系数的含义重点进行了讨论,使用该模型对品种保护申请进行了估计。估计结果显示: 大多数变量系数显著,证明...
  • 第15章Probit回归(概率单位回归).pdf
  • probit_logistic_figure_20210409100939.m
  • 论文研究-网络位置、独立董事治理与公司违规——基于部分可观测Bivariate Probit模型.pdf, 基于社会网络视角研究独立董事治理功能是公司治理研究的最新发展领域. 本文...
  • Probit and Logit Models Stata Program and Output.pdf ,介绍Probit ,logit模型的statad代买,侵权删。
  • glm():family=binomial(link=logit)和(link=probit)、factor、as.factor、relevel

    R语言学习笔记



    R语言中的数据类型之factor因子


    probit回归

    factor()和as.factor()

    这两个没区别。

    #MEPS DATA
    Hexpend<-read.csv("HealthExpend.csv")	#导入数据
    
    # CHECK THE NAMES,DIMENSION IN THE FILE AND LIST THE FRIST
    names(Hexpend)
    dim(Hexpend)
    Hexpend[1:8,]
    attach(Hexpend)
    n<-dim(Hexpend)[1]
    POSEXP<-seq(0,0,length=n)
    for(i in 1:n){
    	if(EXPENDIP[i]!=0)POSEXP[i]=1}
    
    # ALTERNATIVE - FIT A GENERALIZED LINEAR MODEL;
    PosExpglm = glm(POSEXP~GENDER,family=binomial(link=logit))
    summary(PosExpglm)
    logLik(PosExpglm)
    summary(POSEXP)
    
    # FULL LOGIT MODEL
    PosExpglmFull=glm(POSEXP~AGE+GENDER
    	+factor(RACE)+factor(REGION)+factor(EDUC)+factor(PHSTAT)+factor(ANYLIMIT)+factor(INCOME)+factor(insure),
    	family=binomial(link=logit))
    summary(PosExpglmFull)
    logLik(PosExpglmFull)
    
    Gender<-as.factor(GENDER)
    PosExpglmFull=glm(POSEXP~AGE+C(Gender,base=1)
    	+as.factor(RACE)+as.factor(REGION)+as.factor(EDUC)+as.factor(PHSTAT)+as.factor(ANYLIMIT)+as.factor(INCOME)+as.factor(insure),
    	family=binomial(link=logit))
    summary(PosExpglmFull)
    

    relevel()

    # CHANGE REFERANCE LEVELS TO AGREE WITH BOOK(DONE IN SAS)
    RACE=relevel(factor(RACE),ref="WHITE")
    REGION=relevel(factor(REGION),ref="WEST")
    EDUC=relevel(factor(EDUC),ref="LHIGHSC")
    PHSTAT=relevel(factor(PHSTAT),ref="EXCE")
    INCOME=relevel(factor(INCOME),ref="POOR")
    

    这里给的例子里没有factor,但不加会报错。
    可以对比relevel前后的结果

    3
    4

    案例11.4复刻 glm函数

    Regression Modeling with Actuarial and Financial P315 案例11.4 Application: Medical Expenditures

    下载数据HealthExpend.csv

    HealthExpend<-read.csv("HealthExpend.csv")	#导入数据
    str(HealthExpend)
    
    attach(HealthExpend)
    n=2000	#共有2000个数据
    
    sink("result.txt")	#回归结果等储存在这里
    

    整理变量

    ## 整理变量
    #Ethnicity即RACE或RACE1
    ASIAN<-seq(0,0,length=n)	# 1 if Asian 4.3 4.7
    BLACK<-seq(0,0,length=n)	# 1 if Black 14.8 10.5
    NATIVE<-seq(0,0,length=n)	#1 if Native 1.1 13.6???这个的描述统计结果不符
    	#WHITE Reference level 79.9 7.5???other怎么处理
    
    #Region即REGION或REGION1
    NORTHEAST<-seq(0,0,length=n)	# 1 if Northeast 14.3 10.1
    MIDWEST<-seq(0,0,length=n)	# 1 if Midwest 19.7 8.7
    SOUTH<-seq(0,0,length=n)	#1 if South 38.2 8.4
    	#WEST Reference level 27.9 5.4
    
    #Education 即EDUC
    COLLEGE<-seq(0,0,length=n)	# 1 if college or higher degree 27.2 6.8
    HIGHSCHOOL<-seq(0,0,length=n)	# 1 if high school degree 43.3 7.9
    	#Reference level is lower than high school degree 29.5 8.8
    
    #Self-rated physical health 即PHSTAT和PHSTAT1
    POOR<-seq(0,0,length=n)	# 1 if poor 3.8 36.0
    FAIR<-seq(0,0,length=n)	# 1 if fair 9.9 8.1
    GOOD<-seq(0,0,length=n)	# 1 if good 29.9 8.2
    VGOOD<-seq(0,0,length=n)	# 1 if very good 31.1 6.3
    	#Reference level is excellent health 25.4 5.1
    
    
    #Income compared to poverty line 即INCOME和INCOME1
    HINCOME<-seq(0,0,length=n)	# 1 if high income 31.6 5.4
    MINCOME<-seq(0,0,length=n)	# 1 if middle income 29.9 7.0
    LINCOME<-seq(0,0,length=n)	# 1 if low income 15.8 8.3
    NPOOR<-seq(0,0,length=n)	# 1 if near poor 5.8 9.5
    	#Reference level is poor/negative 17.0 13.0
    
    #Insurance coverage insurance in any month of 2003
    INSURE=insure# 1 if covered by public/private health 77.8 9.2 0 if have no health insurance in 2003 22.3 3.1
    
    #因变量
    Y<-seq(0,0,length=n)
    
    for(i in 1:n)
    {
    	if(RACE[i]=="ASIAN")	ASIAN[i]=1
    	if(RACE[i]=="BLACK")	BLACK[i]=1
    	if(RACE[i]=="NATIV")	NATIVE[i]=1
    	
    	if(REGION[i]=="NORTHEAST")	NORTHEAST[i]=1
    	if(REGION[i]=="MIDWEST")	MIDWEST[i]=1	
    	if(REGION[i]=="SOUTH")	SOUTH[i]=1	
    
    	if(EDUC[i]=="COLLEGE")	COLLEGE[i]=1	
    	if(EDUC[i]=="HIGHSCH")	HIGHSCHOOL[i]=1
    
    	if(PHSTAT[i]=="POOR")	POOR[i]=1	
    	if(PHSTAT[i]=="FAIR")	FAIR[i]=1
    	if(PHSTAT[i]=="GOOD")	GOOD[i]=1
    	if(PHSTAT[i]=="VGOO")	VGOOD[i]=1
    	
    	if(INCOME[i]=="HINCOME")	HINCOME[i]=1
    	if(INCOME[i]=="MINCOME")	MINCOME[i]=1
    	if(INCOME[i]=="LINCOME")	LINCOME[i]=1
    	if(INCOME[i]=="NPOOR")	NPOOR[i]=1
    
    	if(EXPENDIP[i]!=0)	Y[i]=1
    }
    
    data<-data.frame(Y,
    	AGE,GENDER,
    	ASIAN,BLACK,NATIVE,
    	NORTHEAST,MIDWEST,SOUTH,
    	COLLEGE,HIGHSCHOOL,
    	POOR,FAIR,GOOD,VGOOD,
    	MNHPOOR,ANYLIMIT,
    	HINCOME,MINCOME,LINCOME,NPOOR,
    	INSURE)	#用于回归的数据集
    detach(HealthExpend)
    attach(data)
    m<-length(data)	#变量数
    summary(data)	#描述统计
    

    回归:Logistic和Probit——glm()

    family=binomial(link=logit)

    Logistic.full<-glm(Y~.,data=data,family=binomial(link=logit))
    print(summary(Logistic.full))
    
    Logistic.reduced<-glm(Y~.-AGE-MNHPOOR,data=data,family=binomial(link=logit))
    print(summary(Logistic.reduced))
    

    family=binomial(link=probit)

    Probit.reduced<-glm(Y~.-AGE-MNHPOOR,data=data,family=binomial(link=probit))
    print(summary(Probit.reduced))
    

    计算log likelihood

    这个可以直接用logLik函数计算,这里是自己手动写了个程序……
    我也是对自己无语了哈哈哈

    #计算log likelihood===============================
    LogLikelihood<-c(0,0,0)
    logit<-function(z){ 
    	1/(1+exp(-z))
    }
    #Logistic.full-----------------------------------
    result<-seq(0,0,length=2000)
    z<-seq(0,0,length=2000)
    pai<-seq(0,0,length=2000)
    
    for(i in 1:n){
    
    	for(j in 1:(m-1)){
    		z[i]<-coef(Logistic.full)[j+1]*data[i,j+1]+z[i]
    		}
    	z[i]=z[i]+coef(Logistic.full)[1]
    	pai[i]=logit(z[i])
    	result[i]=Y[i]*log(pai[i])+(1-Y[i])*log(1-pai[i])
    }
    
    LogLikelihood[1]<-sum(result)
    
    #Logistic.reduced-----------------------------------
    result<-seq(0,0,length=2000)
    z<-seq(0,0,length=2000)
    pai<-seq(0,0,length=2000)
    
    for(i in 1:n){
    
    	for(j in 1:(m-3)){
    		z[i]<-coef(Logistic.reduced)[j+1]*data[i,names(coef(Logistic.reduced)[j+1])]+z[i]
    		}
    	z[i]=z[i]+coef(Logistic.reduced)[1]
    	pai[i]=logit(z[i])
    	result[i]=Y[i]*log(pai[i])+(1-Y[i])*log(1-pai[i])
    }
    
    LogLikelihood[2]<-sum(result)
    
    
    
    #Probit.reduced-----------------------------------
    result<-seq(0,0,length=2000)
    z<-seq(0,0,length=2000)
    pai<-seq(0,0,length=2000)
    
    for(i in 1:n){
    	for(j in 1:(m-3)){
    		z[i]<-coef(Probit.reduced)[j+1]*data[i,names(coef(Probit.reduced)[j+1])]+z[i]
    		}
    	z[i]=z[i]+coef(Probit.reduced)[1]
    	pai[i]=pnorm(z[i])
    	result[i]=Y[i]*log(pai[i])+(1-Y[i])*log(1-pai[i])
    }
    LogLikelihood[3]<-sum(result)
    
    cat("LogLikelihood: ",LogLikelihood)
    
    sink()
    

    结果输出

    #系数
    write.csv(summary(Logistic.full)$coef,file="Logistic_full.csv")
    write.csv(summary(Logistic.reduced)$coef,file="Logistic_reduced.csv")
    write.csv(summary(Probit.reduced)$coef,file="Probit_reduced.csv")
    
    #描述统计
    col.name<-c("min","max","median","mean","sd")
    describe<-matrix(data = NA, nrow = m, ncol = 5, byrow = FALSE,
           dimnames =list(names(data),col.name))
    
    for(i in 1:m){
    	describe[i,1]=min(data[,i])
    	describe[i,2]=max(data[,i])
    	describe[i,3]=median(data[,i])
    	describe[i,4]=mean(data[,i])
    	describe[i,5]=sd(data[,i])
    }
    
    write.csv(data.frame(describe),file="summary.csv")
    

    比较probit和logistic

    compare

    #Probit模型
    curve(exp(x)/(1+exp(x)),xlim=c(-4,4),col="red",ylab="Distribution Function")
    text(2.3,.75,"Logit Case",cex=1.2,col="red")
    arrows(2.3,.79,1.88,.84,length=0.1,angle=20,col="red")
    
    #Logit模型
    curve(pnorm(x),xlim=c(-4,4),col="blue",add=T,ylab="Distribution Function")
    text(0,.95,"Probit Case",cex=1.2,col="blue")
    arrows(-.1,0.9,.6,.8,length=0.1,angle=20,col="blue")
    
    展开全文
  • 我国金融发展对中小企业创新影响的实证研究---基于Probit模型的分析,张强,程莎,中小企业自主创新是经济转型升级的重要驱动因素,如何利用金融支持提升其自主创新能力,为中小企业创新的可持续发展提供保障,是
  • 线性回归、logit回归、probit回归

    千次阅读 2021-04-09 23:28:31
    回归 文章目录 回归 线性回归 古典线性回归模型的假定: OLS的推导与性质 notation 系数求解 标准误 小样本性质 对应检验 拟合优度检验 T 检验 F检验 从似然比角度看F统计量 大样本情况 假定: 性质: Probit&Logistic ...

    回归

    主要借鉴高级计量经济学及Stata应用 第2版_陈强_北京:高等教育出版社_2014.04_669_13526050

    文中所提"书"即是这本中的内容

    线性回归

    Y i = b 0 + b 1 X i + ϵ i , ϵ i ∼ ( 0 , σ 2 ) Y_i=b_0+b_1X_i+\epsilon_i,\qquad \epsilon_i\sim (0,\sigma^2) Yi=b0+b1Xi+ϵi,ϵi(0,σ2)

    Y ^ = b 0 + b 1 X \hat{Y}=b_0+b_1X Y^=b0+b1X

    古典线性回归模型的假定:

    • 总体模型:

      Y i = b 0 + b 1 X i + ϵ i Y_i=b_0+b_1X_i+\epsilon_i Yi=b0+b1Xi+ϵi

    • 严格外生性(strict exogeneity):

      E ( ϵ i ∣ X ) = E ( ϵ i ∣ x 1 , . . . , x n ) = 0 E(\epsilon_i|X)=E(\epsilon_i|x_1,...,x_n)=0 E(ϵiX)=E(ϵix1,...,xn)=0

      即扰动项与所有的解释变量都不相关

    • 不存在严格多重共线性:

      数据矩阵X列满秩 r a n k ( X ) = K rank(X)=K rank(X)=K

    • 球型扰动项

      V a r ( ϵ ∣ X ) = E ( ϵ ϵ ′ ∣ X ) = σ 2 I n Var(\epsilon|X)=E(\epsilon\epsilon'|X)=\sigma^2I_n Var(ϵX)=E(ϵϵX)=σ2In

      条件同方差

    在计量中不满足上述要求会出现以下问题:

    • 自相关性( σ i j ≠ 0 \sigma_{ij}\neq 0 σij=0),即残差之间有相关性,一般存在于时间序列模型中

    • 异方差: σ i ≠ σ j \sigma_i\neq \sigma_j σi=σj

    • 异质性:样本中不同部分不适合共用一个模型

    • 内生性:扰动项与 x i x_i xi相关,对此进行内生性检验

      上面两个问题在大样本中都有一点容忍度,但是内生变量在任何情况都不能存在

    OLS的推导与性质

    OLS(ordinary least square)

    notation
    • SSR(残差平方和) =ESS(explained sum of squares)
    • SSE(回归平方和) =RSS(residual sum of squares)
    • SST(总体偏差平方和)
    系数求解

    to min: S S T = ( y − X β ~ ) T ( y − X β ~ ) SST=(y-X\tilde\beta)^T(y-X\tilde\beta) SST=(yXβ~)T(yXβ~)

    ∂ S S T ∂ β ~ = − 2 ( y T X ) T + ( X T X β ~ ) + ( β ~ T X T X ) T = 0 ⇒ β ~ = ( X T X ) − 1 X T y \frac{\partial SST}{\partial \tilde{\beta}}=-2(y^TX)^T+(X^TX \tilde{\beta})+(\tilde{\beta}^T X^TX)^T=0 \Rightarrow \tilde{\beta}=(X^TX)^{-1}X^Ty β~SST=2(yTX)T+(XTXβ~)+(β~TXTX)T=0β~=(XTX)1XTy

    验证Hessian矩阵: ∂ 2 ( S S R ) ∂ β ~ ∂ β ~ ′ = 2 X T X \frac{\partial^2(SSR)}{\partial\tilde \beta\partial\tilde\beta'}=2X^TX β~β~2(SSR)=2XTX

    PS:可知不存在严格多重共线性保证了 ( X T X ) (X^TX) (XTX)可逆

    标准误

    扰动项方差 σ 2 = V a r ( ϵ i ) \sigma^2=Var(\epsilon_i) σ2=Var(ϵi)

    有无偏估计: s 2 = 1 n − K ∑ i = 1 n e i 2 s^2=\frac{1}{n-K}\sum_{i=1}^ne_i^2 s2=nK1i=1nei2

    因为必须有 X ′ e = 0 X'e=0 Xe=0从而有 n − K n-K nK(样本量-系数个数)个 e i e_i ei相互独立

    s被称为标准误

    小样本性质
    • 线性性 β ~ = b = ( X T X ) − 1 X T y \tilde\beta=b=(X^TX)^{-1}X^Ty β~=b=(XTX)1XTy是y的线性组合

    • 无偏性 E ( b ∣ X ) = β , E ( b ) = β E(b|X)=\beta,E(b)=\beta E(bX)=β,E(b)=β:

      b − β = ( X T X ) − 1 X T ( X β + ϵ ) − β = ( X T X ) − 1 X T ϵ b-\beta=(X^TX)^{-1}X^T(X\beta+\epsilon)-\beta=(X^TX)^{-1}X^T\epsilon bβ=(XTX)1XT(Xβ+ϵ)β=(XTX)1XTϵ

      E ( b − β ∣ X ) = E ( A ϵ ∣ X ) = 0 E(b-\beta|X)=E(A\epsilon|X)=0 E(bβX)=E(AϵX)=0

    • V a r ( b ∣ X ) = σ 2 ( X T X ) − 1 Var(b|X)=\sigma^2(X^TX)^{-1} Var(bX)=σ2(XTX)1

    • Gauss-Markov Theorem:

      OLS是最佳线性无偏估计,即在所有对 β \beta β的线性无偏估计中 V a r ( b ∣ X ) ≤ V a r ( β ^ ∣ X ) Var(b|X)\leq Var(\hat \beta|X) Var(bX)Var(β^X)

      证明具体见书P19

    • E ( s 2 ∣ X ) = σ 2 E(s^2|X)=\sigma^2 E(s2X)=σ2

      见书P20,具体使用了 ϵ T M ϵ = t r ( ϵ T M ϵ ) = t r ( M ϵ ϵ T ) 的 技 巧 \epsilon^TM\epsilon=tr(\epsilon^TM\epsilon)=tr(M\epsilon\epsilon^T)的技巧 ϵTMϵ=tr(ϵTMϵ)=tr(MϵϵT)

    对应检验

    拟合优度检验

    R 2 = ∑ i = 1 n ( y ^ i − y ˉ ) 2 ∑ i = 1 n ( y I − y ^ ) 2 = 1 − S S R S S T ≤ 1 R^2=\frac{\sum_{i=1}^n (\hat y_i-\bar y)^2}{\sum_{i=1}^n(y_I-\hat y)^2}=1-\frac{SSR}{SST}\leq 1 R2=i=1n(yIy^)2i=1n(y^iyˉ)2=1SSTSSR1

    R 2 R^2 R2越高拟合程度越好,其中SST是样本数据固有性质,SSR是拟合效果(增加解释变量必然增大)

    T 检验

    目的:检验系数 β i \beta_i βi是否等于某固定值 β ˉ k \bar \beta_k βˉk

    H 0 : β i = β ^ k ⇔ H 1 : β i ≠ β ^ k H_0:\beta_i=\hat \beta_k\Leftrightarrow H_1:\beta_i\neq \hat \beta_k H0:βi=β^kH1:βi=β^k

    检验统计量: t k = b k − β ˉ k S E ( b K ) = b k − β ˉ k s 2 ( X T X ) k k − 1 ∼ t ( n − K ) t_k=\frac{b_k-\bar \beta_k}{SE(b_K)}=\frac{b_k-\bar \beta_k}{\sqrt{s^2(X^TX)_{kk}^{-1}}}\sim t(n-K) tk=SE(bK)bkβˉk=s2(XTX)kk1 bkβˉkt(nK)

    对于 t k t_k tk服从T分布的证明见书P21,主要用到幂等矩阵M, , X = { x 1 , . . . , x n } x i ∼ N ( 0 , 1 ) , X T M X ∼ r a n k ( M ) ∼ t r a c e ( M ) ,X=\lbrace x_1,...,x_n\rbrace x_i \sim N(0,1), X^T MX\sim rank(M)\sim trace(M) ,X={x1,...,xn}xiN(0,1),XTMXrank(M)trace(M) 相关证明见Appendix

    F检验

    目的:检验模型中m个线性假设是否同时成立,或者有无多余的方程或者自相矛盾的方程。

    在这里插入图片描述

    从似然比角度看F统计量

    在这里插入图片描述

    R 2 R^2 R2代表拟合优度

    大样本情况

    对不用样本残差与解释变量的独立性不做要求,且样本渐进独立且平稳,参数估计渐进一致。

    假定:
    • 线性假定 y i = x i ′ β + ϵ i , ( i = 1 , 2 , . . . , n ) y_i=x_i'\beta+\epsilon_i,\qquad(i=1,2,...,n) yi=xiβ+ϵi,(i=1,2,...,n)
    • { y i , X i } \lbrace y_i,X_i\rbrace {yi,Xi}是渐进独立的平稳过程
    • 前定解释变量: ϵ i \epsilon_i ϵi只需与 X i X_i Xi不相关,和其它样本解释变量可以相关
    • 秩条件: E ( x i x i T ) E(x_ix_i^T) E(xixiT)非退化
    • g i g_i gi为鞅差分序列,且协方差矩阵 S = E ( g i g i ′ ) = E ( ϵ i 2 x i x i ′ ) S=E(g_ig_i')=E(\epsilon_i^2x_ix_i') S=E(gigi)=E(ϵi2xixi)非退化
    性质:

    详见书P58

    • 一致估计 b → p β b\overset{p}\rightarrow \beta bpβ
    • b b b渐进正态
    • A v a r ( b ) Avar(b) Avar(b)有一致估计

    Probit&Logistic

    model

    公式

    P r o b i t : P ( y = 1 ∣ x ) = F ( x , β ) = Φ ( x T β ) Probit:P(y=1|x)=F(x,\beta)=\Phi(x^T\beta) Probit:P(y=1x)=F(x,β)=Φ(xTβ)

    L o g i s t i c : P ( y = 1 ∣ x ) = F ( x , β ) = Λ ( x T β ) = e x p ( x T β ) 1 + e x p ( x T β ) Logistic:P(y=1|x)=F(x,\beta)=\Lambda(x^T\beta)=\frac{exp(x^T\beta)}{1+exp(x^T\beta)} Logistic:P(y=1x)=F(x,β)=Λ(xTβ)=1+exp(xTβ)exp(xTβ)

    对比:Logistic has fat tails ∼ T 7 \sim T_7 T7

    边际效应(marginal effect)

    ∂ P ( y = 1 ∣ x ) ∂ x k \frac{\partial P(y=1|x)}{\partial x_k} xkP(y=1x)

    常用类型

    • 平均边际效应:每个观测值上的边际效应的简单算术平均
    • 样本均值处的边际效应: x = x ˉ x=\bar x x=xˉ处的边际效应
    对logit模型:几率比(odds ratio)

    p = P ( y = 1 ∣ x ) , p 1 − p : 比 率 比 / 相 对 风 险 p=P(y=1|x),\frac{p}{1-p}:比率比/相对风险 p=P(y=1x),1pp:/

    l n ( p 1 − p ) = x T β ln(\frac{p}{1-p})=x^T\beta ln(1pp)=xTβ

    在stata中: e x p ( β ^ i ) exp(\hat\beta_i) exp(β^i)称为第i个变量对应的几率比,表示增加1单位 x i x_i xi, p 1 − p \frac{p}{1-p} 1pp 变化的比值:

    p ∗ 1 − p ∗ p 1 − p = e ( x 1 , . . . , x i + 1 , . . . ) T β e ( x 1 , . . . , x i , . . . ) T β = e β i \frac{\frac{p^*}{1-p^*}}{\frac{p}{1-p}}=\frac{e^{(x_1,...,x_i+1,...)^T\beta}}{e^{(x_1,...,x_i,...)^T\beta}}=e^{{\beta_i}} 1pp1pp=e(x1,...,xi,...)Tβe(x1,...,xi+1,...)Tβ=eβi

    检验

    拟合优度

    借鉴似然函数最大值 R 2 = l n L 0 − l n L 1 l n L 0 R^2=\frac{ln L_0-ln L_1}{lnL_0} R2=lnL0lnL0lnL1

    L 0 : L_0: L0:常数项为唯一解释变量的对数似然函数的最大值

    L 1 : L_1: L1:原模型的对数似然函数的最大值

    Probit:异方差问题

    P ( y i = 1 ∣ x I ) = Φ ( x i T β / σ ) , σ = 1 P(y_i=1|x_I)=\Phi(x_i^T\beta/\sigma),\sigma=1 P(yi=1xI)=Φ(xiTβ/σ),σ=1

    异方差情况下:

    P ( y i = 1 ∣ x I ) = Φ ( x i T β / σ i ) , σ i 2 = e x p ( z i T δ ) P(y_i=1|x_I)=\Phi(x_i^T\beta/\sigma_i),\sigma_i^2=exp(z_i^T\delta) P(yi=1xI)=Φ(xiTβ/σi),σi2=exp(ziTδ)

    对应检验异方差:对 H 0 : δ = 0 H_0:\delta=0 H0:δ=0进行LR检验

    多元回归变量旋转

    依据:F统计量

    在这里插入图片描述

    Algorithm for subset selection

    • Forward selection (逐步加入 p < α p<\alpha p<α的变量)

    • Backward elimination (先对全体进行再逐个删 p > α p>\alpha p>α的变量)

    • Stepwise regression (每次加入后删去 p > α p>\alpha p>α的变量)

    在这里插入图片描述

    Appendix

    幂等矩阵M, , X = { x 1 , . . . , x n } x i ∼ N ( 0 , 1 ) , X T M X ∼ r a n k ( M ) ∼ t r a c e ( M ) ,X=\lbrace x_1,...,x_n\rbrace x_i \sim N(0,1), X^T MX\sim rank(M)\sim trace(M) ,X={x1,...,xn}xiN(0,1),XTMXrank(M)trace(M)

    来自知乎相关问题王赟 Maigo的回答
    在这里插入图片描述

    展开全文
  • 动态 Probit 模型及 Stata 实现

    千次阅读 2019-09-13 09:49:49
    动态 Probit 模型简介 2.1 模型设定 2.2 模型估计 3. Stata应用 3.1 数据简介 3.2 旧方法 (手动处理) 3.3 新方法:外部命令 xtpdyn 3.3.1 下载安装 3.3.2 基本语法 3.2.3 辅助语法 4. 参考文献 关于我们   本篇...

    作者:万源星 (浙江工商大学财务与会计学院)
    邮箱:408887469@qq.com

    连享会 - 与君分享 lianxh.cn


    本篇推文介绍动态面板 Porbit 模型 (Dynamic Probit Model) 的背景、模型及 Stata 应用。

    1. 背景简介

    动态面板数据的一般形式为:

    y i t = α i + ϕ y i t − 1 + β x i t + ε i t y_{it} = \alpha_{i} + \phi y_{it-1} + \beta x_{it} + \varepsilon_{it} yit=αi+ϕyit1+βxit+εit

    模型将 y i t − 1 y_{it-1} yit1 (可以推广到多阶滞后项) 作为解释变量,体现了被解释变量动态变化特征的模型。然而, y i t − 1 y_{it-1} yit1 与随机干扰项 ε i t \varepsilon_{it} εit 具有很强的相关性,这违背了解释变量与随机干扰项不相关的基本假设,即存在内生性问题。基于此,Arellano and Bond (1991) 提出了一阶差分广义矩 (First-difference GMM) 估计法,先通过一阶差分消除个体效应,再用被解释变量的水平值作为工具变量,降低了模型的内生性。Arellano and Bover (1995) 、Blundell and Bond (1998) 在此基础上进一步提出了系统广义矩 (System GMM) ,将被解释变量的差分值作为工具变量,减少了一阶差分广义矩估计法的偏误。

    然而,Probit、Logit 等非线性模型的被解释变量是二元变量 (Binary) 或取值有限的离散变量 (Dichotomous) ,模型用极大似然法进行估计,一阶差分广义矩估计法和系统广义矩估计法难以适用。由此,本文引入动态面板 Porbit 模型 (Dynamic Probit Model) ,下文重点介绍模型的推导过程和 Stata 应用。

    2. 动态 Probit 模型简介

    2.1 模型设定

    假设研究个体行为 (比如上市公司研发) y i t y_{it} yit,即个体 i i i 在第 t t t 年是否执行该行为,它与上一年该行为 y i t − 1 y_{it-1} yit1,以及各种决定因素相关,那么动态面板 Probit 模型可以设定为:

    y i t ∗ = c i + β y i t − 1 + ρ x i t + γ z i t + ε i t (1) y_{it}^{*} = c_{i} + \beta y_{it-1} + \rho x_{it} + \gamma z_{it} + \varepsilon_{it} \tag{1} yit=ci+βyit1+ρxit+γzit+εit(1)

    其中, y i t ∗ y_{it}^{*} yit y i t y_{it} yit 的关系为: y i t ∗ y_{it}^{*} yit 大于临界值 0 时, y i t = 1 y_{it}=1 yit=1,小于等于 0 时, y i t = 0 y_{it}=0 yit=0 c i c_{i} ci 是不可观测的个体效应; x i t x_{it} xit z i t z_{it} zit 为解释变量,前者随时间变化而发生改变,后者则不具备动态变化特征; ε i t \varepsilon_{it} εit 是随机干扰项,这里假设 x i t x_{it} xit z i t z_{it} zit ε i t \varepsilon_{it} εit 不相关,即严格外生性。

    同时,根据 Wooldridge (2005)、Skrondal and Rabe-Hesketh (2014)、Grotti and Cutuli (2018) 的概念及模型,为了消除不可消除的个体效应, c i c_{i} ci 的表达式为:

    c i t = α 0 + α 1 y i 0 + α 2 x ‾ i + α 3 x i 0 + a i (2) c_{it} = \alpha_{0} + \alpha_{1} y_{i0} + \alpha_{2} \overline{x}_{i} + \alpha_{3} x_{i0} + a_{i} \tag{2} cit=α0+α1yi0+α2xi+α3xi0+ai(2)

    其中, y i 0 y_{i0} yi0 x i 0 x_{i0} xi0 是具有动态变化特征变量的初始观测值; x ‾ i \overline{x}_{i} xi 是个体从时间维度上取均值,即 x ‾ i = 1 / T ∑ i = 0 T x i t \overline{x}_{i} = 1/T \sum_{i=0}^{T} x_{it} xi=1/Ti=0Txit a i a_{i} ai 是随机干扰项,服从正态分布。

    2.2 模型估计

    为便于理解,令 δ w i t ∗ = c i + β y i t − 1 + ρ x i t + γ z i t \delta w_{it}^{*} = c_{i} + \beta y_{it-1} + \rho x_{it} + \gamma z_{it} δwit=ci+βyit1+ρxit+γzit,则式 (1) 可改写成:

    y i t ∗ = δ w i t ∗ + ε i t (4) y_{it}^{*} = \delta w_{it}^{*} + \varepsilon_{it} \tag{4} yit=δwit+εit(4)

    假设 F ( ⋅ ) F \left (\cdot \right) F() ε i t \varepsilon_{it} εit 的连续单调递增的概率密度分布函数,即 1 < F ( ⋅ ) < 0 1< F \left (\cdot \right )< 0 1<F()<0,则有:

    P ( y i t = 1 ) = P ( y i t ∗ > 0 ) = P ( ε i t > − δ w i t ∗ ) = 1 − P ( ε i t ≤ − δ w i t ∗ ) = 1 − F ( − δ w i t ∗ ) \begin{aligned} P\left ( y_{it} = 1 \right ) &= P\left ( y_{it}^{*} > 0 \right ) \\ &= P\left ( \varepsilon_{it} > - \delta w_{it}^{*} \right) \\ &=1- P\left ( \varepsilon_{it} \leq - \delta w_{it}^{*} \right) \\ &= 1 - F\left ( -\delta w_{it}^{*} \right) \end{aligned} P(yit=1)=P(yit>0)=P(εit>δwit)=1P(εitδwit)=1F(δwit)

    根据逻辑分布的对称性,整理可得:
    P ( y i t = 1 ) = F ( δ w i t ∗ ) (5) P\left ( y_{it} = 1 \right ) = F\left ( \delta w_{it}^{*} \right) \tag{5} P(yit=1)=F(δwit)(5)

    P ( y i t = 0 ) = 1 − F ( δ w i t ∗ ) (6) P\left ( y_{it} = 0 \right ) = 1 - F\left ( \delta w_{it}^{*} \right) \tag{6} P(yit=0)=1F(δwit)(6)

    根据 (5) 和 (6),可以将式 (1) 改写成:

    y i t = E ( y i t ) + ε i t = p i t + ε i t = 1 − F ( − δ w i t ∗ ) + ε i t (7) y_{it}=E \left ( y_{it} \right ) + \varepsilon_{it} = p_{it} +\varepsilon_{it} = 1 - F\left ( -\delta w_{it}^{*} \right) + \varepsilon_{it} \tag{7} yit=E(yit)+εit=pit+εit=1F(δwit)+εit(7)

    根据逻辑分布的对称性,式 (7) 可简化成:

    y i t = F ( δ w i t ∗ ) + ε i t (8) y_{it} = F\left ( \delta w_{it}^{*} \right) + \varepsilon_{it} \tag{8} yit=F(δwit)+εit(8)

    y i t y_{it} yit 关于它的条件均值的一个回归。如果将 F ( ⋅ ) F \left (\cdot \right ) F() 设定为标准正态分布函数 Φ ( ⋅ ) \Phi \left (\cdot \right ) Φ(),即:

    F ( z ) = Φ ( z ) = ∫ − ∞ z 1 2 π e − t 2 2 d t (9) F \left (z \right )=\Phi \left (z \right ) = \int_{ - \infty }^{ z } {\frac{1}{\sqrt{2\pi }} e^{-\frac{t^{2}}{2}}} dt \tag{9} F(z)=Φ(z)=z2π 1e2t2dt(9)

    由此,式 (5) 可以写成 :

    P ( y i t = 1 ) = F ( δ w i t ∗ ) = Φ ( δ w i t ∗ ) (10) P\left ( y_{it} = 1 \right )=F\left ( \delta w_{it}^{*} \right)= \Phi \left (\delta w_{it}^{*} \right ) \tag{10} P(yit=1)=F(δwit)=Φ(δwit)(10)

    非线性模型的估计方法是极大似然法。那么,假设我们有一组随机样本,为了得到以解释变量为条件的极大似然估计量,可将动态面板 Probit 模型改写成:

    P ( y i t ) = [ F ( δ w i t ∗ ) ] y i t [ 1 − F ( δ w i t ∗ ) ] 1 − y i t , y = 0 , 1 (11) P\left ( y_{it} \right )=\left [ F\left ( \delta w_{it}^{*} \right) \right ]^{y_{it}} \left [ 1- F\left ( \delta w_{it}^{*} \right) \right ]^{1-y_{it}} ,y=0,1 \tag{11} P(yit)=[F(δwit)]yit[1F(δwit)]1yit,y=0,1(11)

    其中, i = 1 , 2.... N , t = 1 , 2.... T i=1,2....N,t=1, 2....T i=1,2....Nt=1,2....T。显然,当 y i t = 1 y_{it}=1 yit=1 时, P ( y i t = 1 ) = F ( δ w i t ∗ ) P\left ( y_{it} = 1 \right )=F\left (\delta w_{it}^{*} \right) P(yit=1)=F(δwit);当 y i t = 0 y_{it}=0 yit=0 时, P ( y i t = 0 ) = 1 − F ( δ w i t ∗ ) P\left ( y_{it} = 0 \right )= 1- F\left ( \delta w_{it}^{*} \right) P(yit=0)=1F(δwit)。因此,根据式 (11) 可获得似然函数:

    L = ∏ i = 1 N ∏ t = 1 T [ F ( δ w i t ∗ ) ] y i t [ 1 − F ( δ w i t ∗ ) ] 1 − y i t (12) L=\prod_{i=1}^{N}\prod_{t=1}^{T}\left [F\left (\delta w_{it}^{*} \right) \right ]^{y_{it}} \left [ 1- F\left (\delta w_{it}^{*} \right) \right ]^{1-y_{it}} \tag{12} L=i=1Nt=1T[F(δwit)]yit[1F(δwit)]1yit(12)

    对数似然函数为:

    l n L = ∑ i = 1 N ∑ t = 1 T { y i t ⋅ l n [ F ( δ w i t ∗ ) ] + ( 1 − y i t ) ⋅ l n [ 1 − F ( δ w i t ∗ ) ] } (13) lnL=\sum_{i=1}^{N} \sum_{t=1}^{T} \left \{ y_{it} \cdot ln\left [ F\left (\delta w_{it}^{*} \right ) \right ]+ (1 - y_{it}) \cdot ln[1-F\left ( \delta w_{it}^{*} \right )]\right \} \tag{13} lnL=i=1Nt=1T{yitln[F(δwit)]+(1yit)ln[1F(δwit)]}(13)

    最大化 l n L lnL lnL 的一阶条件为 ∂ l n L / ∂ δ = 0 {\partial lnL}/{\partial \delta}=0 lnL/δ=0,即:

    ∂ l n L ∂ δ = ∑ i = 1 N ∑ t = 1 T [ δ w i t ∗ ⋅ f ( δ w i t ∗ ) ⋅ y i t − F ( δ w i t ∗ ) F ( 1 − w i t ∗ ) ⋅ F ( δ w i t ∗ ) ] = 0 (14) \frac{\partial lnL}{\partial \delta}=\sum_{i=1}^{N} \sum_{t=1}^{T}\left [ \delta w_{it}^{*} \cdot f\left (\delta w_{it}^{*} \right ) \cdot \frac{y_{it}-F\left ( \delta w_{it}^{*} \right )}{F\left (1- w_{it}^{*} \right )\cdot F\left (\delta w_{it}^{*} \right )}\right ]=0 \tag{14} δlnL=i=1Nt=1T[δwitf(δwit)F(1wit)F(δwit)yitF(δwit)]=0(14)

    其中, F ′ ( δ w i t ∗ ) = f ( δ w i t ∗ ) F^{'}\left (\delta w_{it}^{*} \right )=f\left (\delta w_{it}^{*} \right ) F(δwit)=f(δwit)。由于式 (14) 式是非线性模型,因此需要用迭代法进行求解,这里不再详细展开。

    3. Stata应用

    Stata 官方并未提供直接估计动态 Probit 模型的命令。下面,我们先用手动计算的方法估计该模型,进而使用外部命令 xtpdyn 进行估计。前者便于理解模型背后的原理;后者语法简洁,推荐使用。

    3.1 数据简介

    下文将使用 Grotti and Cutuli (2018) 文中的数据,变量定义如下:

    NameLabel定义
    idIdentification number个体 id
    yearYear of survey调查年份
    poorHousehold poverty家庭条件贫困为 1,否则为 0 (被解释变量)
    ageAge of household head家庭户主年龄
    eduEducation of household head高中以下为 0,高中为 1,以上为 2
    blackRace of household head家族户主是黑人为 1,否则为 0
    empEmployment of household head失业为 0,兼职为 1,全职为 2
    marstatMarital status已婚为 0,单身为 1,其他为 2

    3.2 旧方法 (手动处理)

    旧方法的基本步骤如下:

    • Step1:预先手动计算具有动态变化特征变量的初始值 y i 0 、 x i 0 \color{red}{y_{i0}}、 \color{red}{x_{i0}} yi0xi0,以及个体均值 x ‾ i \color{red}{\overline{x}_{i}} xi,去除不可观测的个体效应;
    • Step2:使用 xtprobit}meprobit 命令对动态数据模型进行估计;
    • 以 Grotti and Cutuli (2018,[PDF]) 的例子,命令如下:
    ******************************************* 
    *      手动估计动态面板Probit模型 
    * By-hand replication of model estimation            
    *******************************************
    
    *-Step1:手动计算相关初始值、个体均值
    
      *-数据概况
        use poverty.dta,clear
        xtset id year
        xtdes
      *-去掉缺失值
        gen sample = !mi(id, year, poor, black, age, edu, emp, marstat)
      *-计算具有动态变化特征变量的初始值
        local i "poor age emp marstat"
        foreach var of varlist `i' {
        bys sample id (year): gen `var'__0 = `var'[1] ///
        if sample & sample[1]==1	
      	}
      *-计算具有动态变化特征变量的个体均值
        bys sample id : egen m__age = mean(age)
        bys sample id : egen m1__emp = mean(emp==1)
        bys sample id : egen m2__emp = mean(emp==2)
        bys sample id : egen m1__marstat = mean(marstat==1)
        bys sample id : egen m2__marstat = mean(marstat==2)
    
    *-Step2:估计模型 -xtprobit- -meprobit-
    	
      *-用 xtprobit 估计动态数据模型
        xtset id year
        xtprobit poor iL.poor black c.age i.edu i.emp i.marstat ///
      		i.poor__0 age__0 i.emp__0 i.marstat__0 ///
      		m__age m1__emp m2__emp m1__marstat m2__marstat
        est store xtp
      	
      *-用 meprobit 估计动态数据模型 
        meprobit poor iL.poor black c.age i.edu i.emp i.marstat ///
      		i.poor__0 age__0 i.emp__0 i.marstat__0 ///
      		m__age m1__emp m2__emp m1__marstat m2__marstat || ///
    		id : , intpoints(12)
        est store mep
    

    3.3 新方法:外部命令 xtpdyn

    Grotti and Cutuli (2018, Stata Journal,[PDF]) 详细介绍了 xtpdyn 命令的原理及用法。

    3.3.1 下载安装

    在 Stata 命令窗口输入如下命令,即可自动安装 xtpdyn 命令:

    ssc install xtpdyn, replace
    

    3.3.2 基本语法

    xtpdyn 可以用简单的命令生成变量的初始值及个体均值,语法如下:

    xtpdyn depvar indepvars [if] [in] [weight], uh(varlist)  [re(re_options)]
    
    • depaver 是被解释变量,即模型 (1) 中的 y i t ∗ y_{it}^{*} yit
    • indepvars 是解释变量,即模型 (1) 中的 x i t 和 z i t x_{it} 和 z_{it} xitzit
    • uh() 中放入具有动态变化特征的解释变量,即模型 (1) 中的 x i t x_{it} xit;
    • re() 是指随机效应,如果考虑异方差,可写成 re(vce(robust))

    一个例子 :下面的例子中,xtpdyn 列示了解释变量,以及 poor__0age__0 等变量初始值,m1__empm2__emp 等个体均值的回归结果。

    net get xtpdyn.pkg   // 下载数据
    use "poverty.dta", clear 
    xtset id year
    xtpdyn poor black c.age i.edu i.emp i.marstat, ///
                uh(i.emp i.marstat age) re(vce(robust))
    
    Iteration 0:   log pseudolikelihood = -2101.8697  (not concave)
    …… (output omitted) ……
    Iteration 7:   log pseudolikelihood =  -1963.425  
    
    Mixed-effects probit regression                 Number of obs     =      6,173
    Group variable:              id                 Number of groups  =      1,112
    
                                                    Obs per group:
                                                                  min =          3
                                                                  avg =        5.6
                                                                  max =          9
    
    Integration method: mvaghermite                 Integration pts.  =         12
    
                                                    Wald chi2(20)     =    1154.68
    Log pseudolikelihood =  -1963.425               Prob > chi2       =     0.0000
                                       (Std. Err. adjusted for 1,112 clusters in id)
    --------------------------------------------------------------------------------
                   |               Robust
              poor |      Coef.   Std. Err.      z    P>|z|     [95% Conf. Interval]
    ---------------+----------------------------------------------------------------
                   |
            L.poor |
             Poor  |   .5744007   .0813891     7.06   0.000      .414881    .7339205
                   |
             black |   .5225842   .0688626     7.59   0.000     .3876159    .6575525
               age |  -.0649581   .0129398    -5.02   0.000    -.0903197   -.0395966
                   |
               edu |
      High School  |  -.2264437   .0768907    -2.95   0.003    -.3771466   -.0757407
     More than HS  |  -.5724683   .0960026    -5.96   0.000    -.7606299   -.3843067
                   |
               emp |
        Part-time  |  -.4357227   .1355535    -3.21   0.001    -.7014027   -.1700426
        Full-time  |   -1.23497   .1408088    -8.77   0.000     -1.51095   -.9589893
                   |
           marstat |
           Single  |  -.6937395   .2385683    -2.91   0.004    -1.161325   -.2261542
      Div/sep/wid  |  -.5112024   .1342855    -3.81   0.000    -.7743972   -.2480076
                   |
         1.poor__0 |    .923642    .095303     9.69   0.000     .7368516    1.110432
                   |
            emp__0 |
                1  |   -.148413   .1406438    -1.06   0.291    -.4240698    .1272438
                2  |  -.1200204   .1545567    -0.78   0.437    -.4229458    .1829051
                   |
        marstat__0 |
                1  |  -.3251233   .1983599    -1.64   0.101    -.7139016     .063655
                2  |  -.2616084   .1915875    -1.37   0.172     -.637113    .1138962
                   |
            age__0 |  -.0311116   .0313137    -0.99   0.320    -.0924852    .0302621
           m1__emp |   .2515403   .2492729     1.01   0.313    -.2370256    .7401062
           m2__emp |  -.1801458    .261989    -0.69   0.492    -.6936349    .3333432
       m1__marstat |     1.0849   .3453074     3.14   0.002     .4081096     1.76169
       m2__marstat |   1.051389   .2589864     4.06   0.000     .5437854    1.558993
            m__age |    .081758   .0334461     2.44   0.015     .0162049     .147311
             _cons |  -.1005763    .244535    -0.41   0.681    -.5798561    .3787036
    ---------------+----------------------------------------------------------------
    id             |
         var(_cons)|   .3355826   .0634519                      .2316618     .486121
    --------------------------------------------------------------------------------
    

    3.2.3 辅助语法

    第一xtpdyn 可一键删除或保留新增的变量

     xtpdyn [, keep|drop]
    

    其中,

    • keep 保留变量初始值、个体均值等新增的变量;
    • drop 删除已保留的变量初始值、个体均值等新增的变量;

    第二probat 命令可用于显示动态面板 Probit 模型估计后的事件发生概率 (Transition Probabilities) 及其他统计信息。

     probat, {stats[(atspec)]|prdistr} 
             [margins(margins_options) nq(#)
              showfreq plot keep]
    
    • stats() 选项

      • 记录进入概率 P r ( 1 ∣ 0 ) Pr(1|0) Pr(10)、退出概率 P r ( 0 ∣ 1 ) Pr(0|1) Pr(01)、持续概率 P r ( 1 ∣ 1 ) Pr(1|1) Pr(11)、稳态概率 [ P r ( 1 ∣ 0 ) / P r ( 1 ∣ 0 ) + P r ( 0 ∣ 1 ) ] \left [ Pr\left ( 1|0 \right )/ Pr\left ( 1|0 \right ) +Pr\left ( 0|1 \right )\right ] [Pr(10)/Pr(10)+Pr(01)]、稳态持久度 [ 1 / P r ( 0 ∣ 1 ) ] [1/Pr(0|1)] [1/Pr(01)] 等统计信息;
      • 输入命令 probat, stats 将提供整个数据的统计结果;
      • 输入命令 probat, stats\left (varname1= n \right ) 将提供变量 varname1 的统计信息;
      • 不能与 prdistr 一起使用;
    • prdistr 选项

      • 计算不同水平 c i c_{i} ci 的预期概率 P r ( 1 ∣ 0 ) Pr(1|0) Pr(10) P r ( 1 ∣ 1 ) Pr(1|1) Pr(11) 等。 c i c_{i} ci 产生差异的原因在于,具有个体间的 y i 0 y_{i0} yi0 x ‾ i \overline{x}_{i} xi x i 0 x_{i0} xi0 数值不同;
      • 通过两个维度区分 c i c_{i} ci 的不同水平:第一个维度是被解释变量的初始观测值 y i 0 y_{i0} yi0;第二个维度是根据 x ‾ i \overline{x}_{i} xi x i 0 x_{i0} xi0 的样本分布 (分位数概念)
    • magins 选项: 计算模型的边际效应

    • nq() 选项: 设置计算预期概率时的分位数水平,仅允许 2、3、4、5 和 10 分位数。系统默认为 5 分位数。

    • showfreq 选项:显示不同水平 c i c_{i} ci 的样本频数

    • plot 选项:图示不同水平 c i c_{i} ci 的预期概率结果

    Stata 范例 1:

    • 命令含义是在 3 分位数水平上列示样本分布及 prdistr 的统计结果;
    • 第一张表列示样本分布;
    • 第二张表列示统计结果,以第一行为例:在滞后一期的家庭条件为非贫穷 ( L . p o o r = 0 ) \left ( L.poor=0 \right ) (L.poor=0),家庭条件最初为非贫穷 ( p o o r 0 = 0 ) \left ( poor_0 = 0 \right ) (poor0=0),年龄、雇佣方式等联合均值在第一分位数上 ( u h q = 1 ) \left ( uh_q=1 \right) (uhq=1) 的情况下,家庭条件由贫穷变为非贫穷的概率是 5.26%。
    probat, prdistr nq(3) showfreq 
    
    -> poor__0 = 0
    --------------------------------------------------------------------------
              |            3 quantiles of uhi and Household poverty           
    poor at   | -------- 1 -------    -------- 2 -------    -------- 3 -------
    time t-1  | Not poor      Poor    Not poor      Poor    Not poor      Poor
    ----------+---------------------------------------------------------------
            0 |    1,595       104       1,077        95       1,131        94
            1 |       97        55          81        89          73        87
    --------------------------------------------------------------------------
    -> poor__0 = 1
    --------------------------------------------------------------------------
              |            3 quantiles of uhi and Household poverty           
    poor at   | -------- 1 -------    -------- 2 -------    -------- 3 -------
    time t-1  | Not poor      Poor    Not poor      Poor    Not poor      Poor
    ----------+---------------------------------------------------------------
            0 |       71        21         152        61         115        57
            1 |       53        74         142       351         132       366
    --------------------------------------------------------------------------
    
    --------------------------------------------------------------------------------
    L.poor#poor__0#uh_q |    Prob.     Std. Err.     P>|z|     Lower CI    Upper CI 
    --------------------+-----------------------------------------------------------
                  0 0 1 |    0.0526    0.005130       0.000      0.0425      0.0626 
                  0 0 2 |    0.0815    0.006289       0.000      0.0692      0.0939 
                  0 0 3 |    0.0972    0.008282       0.000      0.0809      0.1134 
                  0 1 1 |    0.2935    0.027079       0.000      0.2404      0.3466 
                  0 1 2 |    0.4713    0.025980       0.000      0.4203      0.5222 
                  0 1 3 |    0.5119    0.028153       0.000      0.4568      0.5671 
                  1 0 1 |    0.1221    0.017207       0.000      0.0884      0.1559 
                  1 0 2 |    0.1707    0.019060       0.000      0.1334      0.2081 
                  1 0 3 |    0.1950    0.021376       0.000      0.1531      0.2369 
                  1 1 1 |    0.4743    0.028741       0.000      0.4180      0.5307 
                  1 1 2 |    0.6440    0.018173       0.000      0.6083      0.6796 
                  1 1 3 |    0.6843    0.018109       0.000      0.6488      0.7198 
    --------------------------------------------------------------------------------
    

    Stata 范例 2:

    • 图为例子一命令中增加 plot 后的结果;
    • 可读信息一:在控制不可观测的个体效应的基础上,家庭条件初始贫穷时,本期继续贫穷的概率更大;
    • 可读信息二:在控制不可观测的个体效应的基础上,上一期家庭条件贫穷时,本期继续贫穷的概率更大;

    4. 参考文献

    • Grotti, R. and Cutuli, G. Xtpdyn: A community-contributed command for fitting dynamic random-effects probit models with unobserved heterogeneity. The Stata Journal, 2018, 18(4): 844-862. [PDF]
    • Rabe-Hesketh, S. and Skrondal, A. Avoiding biased versions of Wooldridge’s simple solution to the initial conditions problem. Economics Letters, 2013, 120: 346-349. [PDF]
    • Skrondal, A. and S. Rabe-Hesketh, S. Handling initial conditions and endogenous covariates in dynamic/transition models for binary data with unobserved heterogeneity. Journal of the Royal Statistical Society, 2014, 63©: 211–237. [PDF]
    • Wooldridge, J. M. Simple solutions to the initial conditions problem in dynamic, nonlinear panel data models with unobserved heterogeneity. Journal of Applied Econometrics, 2005, 20: 39-54. [PDF]


    关于我们

    • Stata连享会 由中山大学连玉君老师团队创办,定期分享实证分析经验。直播间 有很多视频课程,可以随时观看。
    • 你的颈椎还好吗? 您将 ::连享会-主页::::连享会-知乎专栏:: 收藏起来,以便随时在电脑上查看往期推文。
    • 公众号推文分类: 计量专题 | 分类推文 | 资源工具。推文分成 内生性 | 空间计量 | 时序面板 | 结果输出 | 交乘调节 五类,主流方法介绍一目了然:DID, RDD, IV, GMM, FE, Probit 等。
    • 公众号关键词搜索/回复 功能已经上线。大家可以在公众号左下角点击键盘图标,输入简要关键词,以便快速呈现历史推文,获取工具软件和数据下载。常见关键词:
      • 课程, 直播, 视频, 客服, 模型设定, 研究设计,
      • stata, plus,Profile, 手册, SJ, 外部命令, profile, mata, 绘图, 编程, 数据, 可视化
      • DID,RDD, PSM,IV,DID, DDD, 合成控制法,内生性, 事件研究
      • 交乘, 平方项, 缺失值, 离群值, 缩尾, R2, 乱码, 结果
      • Probit, Logit, tobit, MLE, GMM, DEA, Bootstrap, bs, MC, TFP
      • 面板, 直击面板数据, 动态面板, VAR, 生存分析, 分位数
      • 空间, 空间计量, 连老师, 直播, 爬虫, 文本, 正则, python
      • Markdown, Markdown幻灯片, marp, 工具, 软件, Sai2, gInk, Annotator, 手写批注
      • 盈余管理, 特斯拉, 甲壳虫, 论文重现
      • 易懂教程, 码云, 教程, 知乎

    展开全文
  • 为完整掌握电力供应网络中用户节点的跨越调度行为,设计基于WOE-Probit逐步回归的用户跨域行为模式挖掘系统。按照数据挖掘框架的结构类型,连接用户行为处理单元与信息存储模块,实现挖掘系统的硬件执行环境设计。...
  • 请问Python能做probit回归分析吗?Python能做类似多元回归的计量模型吗?自己假设的。
  • spss数据分析常用数据集:probit_LD50.sav 统计分析及模型构建中常用的数据集; 学习软件的时候,会苦于没有数据进行实操,而其实一般分析软件都会自带数据,现在介绍如何获取SPSS软件自带的数据。 纽约时报的一篇...
  • probit回归:即概率单位回归,主要用来测试分析刺激强度与反应比例之间的关系,例如对于指定数量的病人,分析他们的给药剂量与治愈比例之间的关系,此方法运用的典型例子是分析杀虫剂浓度和杀死害虫数量之间的关系,...
  • Probit回归模型.pdf

    2021-12-20 01:01:32
    Probit回归模型.pdf
  • X = wine.iloc[:, :-1] Y = wine['quality'] binary_Y = [] for i in range(len(Y)): if Y[i] (0) else: binary_Y.append(1) Probit probit_model = Probit(binary_Y, sm.add_constant(X)) result = probit_model....
  • 针对Probit分位回归在参数随机化条件下的建模问题,提出基于Metropolis- Hastings算法的贝叶斯Probit分位回归模型.通过分析Probit分位回归模型结构,选择模型的先验分布,运用M-H算法进行参数估计.利用Monte Carlo仿真...
  • 笔记:二元Probit与Logit模型

    千次阅读 2020-01-14 11:00:42
    选择的概率分布常用的是标准正态分布和逻辑分布,相应地形成了两种最常用的二元选择模型——Probit模型与Logit模型。 这两种分布都是对称的,所以 P ( Y i = 1 ) = P ( Y i ∗ > 0 ) = P ( μ i ∗ > − X i β )...
  • Probit和Logit回归

    2021-03-25 20:16:38
    https://blog.csdn.net/weixin_41677876/article/details/106686250
  • 本篇再介绍一种常见的广义线性模型:Logistic模型。该模型主要针对分类结果进行建模。与之功能类似的另一个模型是Probit模型,但较少应用。Logistic模型的形式两点分布,又称伯努...
  • Logit模型和Probit模型

    千次阅读 2021-08-03 21:28:00
    导读 原理 总结 复现论文
  • 如果大家想要进一步深究 Logit 与 Probit 模型的异同,可以参考 Stata 连享会往期推文「 二元选择模型:Probit 还是 Logit?」。 2. 问题与方法 当我们运用 Logit 回归的时候,一般会有两种情况造成回归结果出现以下...
  • 用于 I-II 期剂量寻找试验的贝叶斯双变量 Probit 模型 传统上,在剂量寻找试验中,研究人员采用基于规则的方法进行单独的 1 期和 2 期试验,以寻找最大耐受剂量水平并评估实验药物的功效。 或者,可以利用基于单一...

空空如也

空空如也

1 2 3 4 5 ... 20
收藏数 1,440
精华内容 576
关键字:

probit

友情链接: pintu.rar