1、一、根据美国人口从1790年到1990年间的人口数据(如下表),确定人口指数增长模型和Logistic模型中的待定参数,估计出美国2010年的人口,同时画出拟合效果的图形。表1 美国人口统计数据年 份1790180018101820183018401850人口(106)3.95.37.29.612.917.123.2年 份1860187018801890190019101920人口(106)31.438.650.262.976.092.0106.5年 份193019401950196019701980人口(106)123.2131.7150.7179.3204.0226.5提示:指数增长模型:
2、解:设为第t年美国人口数,记时刻(t=0)的人口为。假设人口增长率为常数r,则有令,有关于的微分方程 ,故可解出得到热口指数增长模型:,(1)将(1)式去对数,可得用MATLAB软件计算可得r=0.214/10年,=-36.6198,预测估计美国2010年的人口数为x(2010)= 598.4,拟合图用MATLAB解题具体步骤如下:程序:t=1790:10:1980;y=3.9 5.3 7.2 9.6 12.9 17.1 23.2 31.4 38.6 50.2 62.9 76 92 106.5 123.2 131.7 150.7 179.3 204 226.5v=log(y);p=polyfi
3、t(t,v,1) %得到指数模型参数值r和a结果:y = Columns 1 through 9 3.9000 5.3000 7.2000 9.6000 12.9000 17.1000 23.2000 31.4000 38.6000 Columns 10 through 18 50.2000 62.9000 76.0000 92.0000 106.5000 123.2000 131.7000 150.7000 179.3000 Columns 19 through 20 204.0000 226.5000p = 0.0214 -36.6198程序:ti=1790:2:1980;vi=polyv
4、al(p,ti);yi=exp(vi)plot(t,y,*,ti,yi,r)legend(数据点(x,y),拟合曲线y=f(x) x(2010)=exp(0.0214*2010-36.6198) %估计2010年美国人口Logistic模型:首先我们可用微积分将连续的微分方程离散化,不妨设x(t)表示第t次所得人口数,r为人口固有增长率,根据logistic模型,为人口容量,以下用k代替,可得: 这里可以假设,则现在我们可以用线性拟合,借助matlab来进行运算得到r,ky=3.9,5.3,7.2,9.6,12.9,17.1,23.2,31.4,38.6,50.2,62.9,76.0,92.0
5、,106.5,123.2,131.7,150.7,179.3,204.0,226.5;t=1:20;for k=1:19 x(k)=(y(k+1)-y(k)/y(k);endA=polyfit(y(1:19),x,1) %求出参数r与sA = -0.0014 0.3368这里求出r=0.3368,k=0.3368/0.0014=240.57143(下面的b0就是这里求出)function y=lin(b,t);f(x)=b(1)./1+(b(1)./3.9-1).*exp(-b(2).*t);t=1:20;y=3.9,5.3,7.2,9.6,12.9,17.1,23.2,31.4,38.6,5
6、0.2,62.9,76.0,92.0,106.5,123.2,131.7,150.7,179.3,204.0,226.5;b0=240.57143,0.3368;b=nlinfit(t,y,lin,b0)b0 = 240.5714 0.3368b = 312.6056 0.2615这里就求出了函数的值r=0.2615,=312.6065,=3.9(最初的值,1790年的)则可求出t=1790:10:1980;y=3.9,5.3,7.2,9.6,12.9,17.1,23.2,31.4,38.6,50.2,62.9,76.0,92.0,106.5,123.2,131.7,150.7,179.3,2
7、04.0,226.5;x1=1:20;F=312.6056./1+80.155282.*exp(-0.2615.*x1);plot(t,y,r*);hold onplot(x1,F,b-);hold off,legend(数据点(xi,yi),拟合曲线y=f(x)xlabel(t),ylabel(x(t)下面是图形预测2010年,当t=23时,求出美国2010年的人口为261.4142结论:从以上两个图形可以看出,logistic拟合的更好,误差更少,预测用logistic模型更加准确。二、f(x)的定义如下:1、写一个函数文件f(x)实现该函数,要求参数x可以是向量;程序:function
8、f=f(x)if x0&x=-4 f=x.2+x-6;elseif x=0&x=2&x=3 f=x.2-5*x+6;else f=x.2-x-1end2、作出该函数的图形;程序:x=-10:0.1:0;y1=x.2+x-6;i=find(x=-4);y1(i)=NaN;plot(x,y1,r)hold onx=0:0.1:10;y2=x.2-5*x+6;i=find(x=2&x=3);y2(i)=NaN;plot(x,y2,g)hold onx=10:0.1:20;y3=x.2-x-1;plot(x,y3,b);legend(y1=x.2+x-6,y2=x.2-5*x+6,y3=x.2-x-1
9、)3、求出f(x)的零点与最值。程序:x=fminbnd(f,-10,20)结果:x = -0.5000即f(x)的零点x=-0.5,最小值为f(x)=-6.25,最大值为无穷。三、财政收入预测问题:财政收入与国民收入、工业总产值、农业总产值、总人口、就业人口、固定资产投资等因素有关。下表列出了1952-1981年的原始数据。1、编程将数据读入MATLAB,对数据进行探索性分析; 2、请你分析表中是否存在不合理数据,如果存在,找出不合理数据; 3.、构造预测模型;4、将以上分析组织成一篇数学建模论文。 表3 1952-1981数据年份国民收入(亿元)工业总产值(亿元)农业总产值(亿元)总人口(
10、万人)就业人口(万人)固定资产投资(亿元)财政收入(亿元)1952598349461574822072944184195358645547558796213648921619547075204916026621832972481955737558529614652232898254195682571555662828230181502681957837798575646532371113928619581028123559865994266002563571959111416815096720726173338444196010791870444662072588038050619617571
11、156434658592559013827119626779644616729525110662301963779104651469172266408526619649431250584704992773612932319651152158163272538286701753931966132219116877454229805212466196712491647697763683081415635219681187156568078534319151273031969137221016888067133225207447197016382747767829923443231256419711
12、780315679085229356203556381972183333657898717735854354658197319783684855892113665237469119741993369689190859373693936551975212142549329242138168462692197620524309955937173883444365719772189492597194974393774547231978247555901058962593985655092219792702606511509754240581564890198027916592119498705418
13、96568826198129276862127310007273280496810财政收入分析与预测问题【摘要】国家的财政收入与国民收入、工业总产值、农业总产值、总人口、就业人口、固定资产投资等因素有关。首先,我们根据所给数据,对数据进行描述性分析。之后,我们对数据进行了回归分析,构造了多元线性回归模型,并获得了模型的回归系数估计值及其置信区间。然后,考虑到每个回归系数置信区间包含零点与否的情况和变量之间存在多重共线性,我们用逐步回归法剔除了工业总产值、就业人口因素变量对模型进行了改进。又考虑到数据的时间序列属性,我们对模型进行了自相关性诊断,作出残差散点图,初步判定其大部分点落在1,3象限,
14、随机误差表现出正自相关趋势。但在之后的D-W检验中,我们计算出了DW值,自相关系数估计值,依照样本容量和回归变量数目,查阅了D-W分布表,得到检验的临界值dL和dU。在分析DW所在区间时,我们发现模型的自相关状态不能确定。但我们将1952年-1980年的各项经济指标数据代入模型,得出的预测值与实际值相当吻合。模型拟合效果非常好最后,利用该模型对1990年和2000年的财政收入作出预测,并对结果进行了分析。关键词:MATLAB 多元线性回归模型 回归分析 逐步回归法一、问题重述国家的财政收入与国民收入、工业总产值、农业总产值、总人口、就业人口、固定资产投资等因素有关,根据题目所给数据,对数据进行
15、分析,剔除不合理数据,构造多元线性回归模型,并利用该模型对1990年和2000年的财政收入作出预测。二、问题假设1、 题目所给的数据准确无误;2、 财政收入为研究指标Y,国民收入、工业总产值、农业总产值、总人口、就业人口、固定资产投资分别自变量;3、 Y与各自变量呈线性函数关系;4、 Y是随机变量,服从均值为零是正态分布,即 三、符号说明:财政收入:国民收入:工业总产值:农业总产值:总人口:就业人口:固定资产投资:与想对应的回归系数, 四、问题分析与模型建立1、问题分析 首先将数据作初步分析。分别将数据导入MATLAB软件中并作出财政收入与6个因素的回归拟合图,并且求出各自拟合优度,F值,p值
16、。如下图(1)到图(7)所示注:MATLAB软件中F(1,29)的临界值为0.004,F(3,25)=0.11582图(1)通过MATLAB软件计算得出=0.9512 F=0.4455 p=0.000表示财政收入有95.12%被国民收入所解释,且P小于0.0001则方程回归效果好,国民收入与财政收入有显著的线性相关。图(2) 通过MATLAB软件计算得出=0.9231 F=0.336 p=0.000表示财政收入有92.31%被国民收入所解释,且P小于0.0001则方程回归效果好,工业总产值与财政收入有显著的线性相关。图(3)通过MATLAB软件计算得出=0.8370 F=0.1438 p=0.
17、000表示财政收入有83.7%被国民收入所解释,且P小于0.0001则方程回归效果好,农业总产值与财政收入有显著的线性相关。图(4)通过MATLAB软件计算得出=0.8756 F=0.1971 p=0.000表示财政收入有87.56%被国民收入所解释,且P小于0.0001则方程回归效果好,总人口与财政收入有显著的线性相关。图(5) 由图(5)可知,1981年就业人口数据存在不合理性,因此将1981年就业人口数据和1981年其他相应的数据剔除,再将剩余的数据作出就业人口与财政收入的拟合图,如图(5-1) 图(5-1) 通过MATLAB软件计算得出=0.8160 F=0.0045 p=0.000表
18、示财政收入有81.6%被国民收入所解释,且P小于0.0001则方程回归效果好,国民收入与财政收入有显著的线性相关。图(6)从上所有拟合图可知,财政收入分别与国民收入、工业总产值等6要素存在着高度线性相关关系,因此,可建立线性回归模型。2.多元线性回归模型基本概念(1)多元线性模型变量Y与其他有关变量的关系称为“回归模型”,其中为均值为0的随机变量。模型中的变量Y称为因变量,通常建立回归模型的目的是要说明因变量的变化规律,并对其进行预测;模型中的称为自变量或“解释变量”,用来说明或解释因变量。当f为线性函数时,回归模型:称为线性回归模型,特别的,当k=1时称为一元线性回归模型。是线性回归的一般形
19、式,其中是未知的参数;是不可观测的随机变量,称为误差项,假定。如果有n次独立的观测数据 ,则线性回归模型可表示成如下形成: 式中,相互独立且服从N(0,)分布。简写成矩阵形式:式中,(2)多元线性回归模型的基本假定假设1:回归模型是正确设定的;假设2:解释变量是非随机的,且各之间不存在严格线性相关性(无完全多重共线性);假设3:各解释变量再所抽取的样本中具有变异性,且随着样本点容量的无限增加,各解释变量的样本方差趋于一个非零常数,即时假设4:随机误差限具有条件零均值、同方差及不序列相关性;假设5:解释变量与随机项不相关;假设6:随机项满足正态分布。(3)参数与的估计 记的估计为,那么我们就可以
20、得到一个k元线性方程,对每一个样本点,可求出相应的值,则为预测值,要使一切的与之间的偏差达到最小,根据最小二乘法的原理,即要求所以只要使达到极小的。由于Q是的一个非负二次型,故其极小值必存在,根据微积分的理论知道,只要Q对的一阶偏导数为0:经整理即得关于的一个线性方程组: (1)其解称为的最小二乘估计。若记(1)的系数矩阵为A,常数项矩阵为B,则A恰为,B恰为: 因而(1)用矩阵形式即为:,在回归分析中存在,这是最小二乘估计,当我们求得了的最小二乘估计后,就可以建立回归方程,从而可以利用它对指标进行预报和控制。为求的估计,引入几个名词。称实测值与回归值的差-为残差,预测值与实测值的关系是:式中
21、,。称为残差向量,其中I为n阶单位矩阵。而称为残差平方和。若rank(X)=k+1n,均方残差(MSE):即为误差的方差的无偏估计。(4)利用回归模型进行预测利用回归方程对因变量的取值进行点预测,点预测是根据回归方程代入自变量的值,得到对应因变量的预测值。回归方程为当自变量的一组新的观测值对应的因变量的预测值为: 3.模型建立根据模型中影响因素的多少及影响因素与预测对象之间的相互关系,回归分析可分为一元线性回归分析、多元线性回归分析和非线性回归分析。本文采用多元线性回归分析。回归分析方法即根据相互影响、相互关联的两个或多个因素(又称为变量)的数据,由不确定的函数关系建立数学模型,确定参数,从而
22、建立函数关系的过程。而将回归分析用于预测和研究未来,寻找待测对象与影响因素之间的数学关系,并采用数学模型予以表达,然后通过确定未来影响因素,间接导出待测数据的过程,称为回归分析预测方法。多元线性回归模型中系数的求解或估计是多元线性回归分折中的基本问题,采用统计学中最成熟使用最广泛的最小二乘估计方法,应用MATLAB软件实现多元线性回归分析,得到多元线性预测模型。从以上的各个一元回归图及对初步的回归分析,我们再引入一个常量回归系数,作出了初步的模型: (1)其中,,是7个待估计的参数,称为总体回归参数;是随机变量,其均值为0,方差为,即服从同一正态分布。4、模型的求解 首先我们剔除掉因为1981
23、年就业人口对财政收入影响异常的特殊点(见图5),之后利用MATLAB统计工具箱中命令regress求解,得到模型(1)的回归系数估计值及其置信区间(置信水平=0.05)、检验统计量R2,F,p的结果见表1。参数参数估计值参数置信区间0-15.5344-366.5816 335.512710.51000.2301 0.78982-0.0259-0.0769 0.02513-0.5905-0.9901 -0.190840.0113-0.0028 0.02545-0.0230-0.0492 0.003260.3419-0.0387 0.7225R2=0.9840,F=225.8953,p=0.000
24、0表1 模型(1)的计算结果结果分析:由表1显示,模型的拟合优度值R2=0.9840表示因变量y(财政收入)的98.40%可由模型确定,模型的整体显著水平F值远远超过F检验的临界值,模型的整体拟合相伴概率=0远小于,因而模型(1)从整体来看是可用的。表1的回归系数给出了模型(1)中的估计值,即,。检查它们的置信区间发现,的置信区间包含零点,但无法确定其6个因素之间的关系,因此不能简单地剔除没通过显著性检验的因素。5、模型改进:5.1改进模型的多重共线性由于计量经济学中多元线性回归模型的经典假设3已经假定变量之间不存在严格相关性,而本文采用时间序列数据作样本,以简单线性形式建立的计量经济学模型,
25、往往存在多重共线性。因此,为排除变量间存在多重共线性,用MATLAB软件进行逐步回归剔除一些变量,从而达到排除变量之间存在多重共线性问题,可得到回归最优模型,结果如图(7) 图(7)注:MATLAB软件中逐步回归法得到的图表中红色点表示剔除的变量,因而可以变量间的排除多重共线性,得到最优模型 从图中数据可知,只有系数通过t检验,利用MATLAB统计工具箱中命令regress求解,得到模型(1)的回归系数估计值及其置信区间(置信水平=0.05)、检验统计量R2,F,p的结果见表2。参数参数估计值参数置信区间163.101079.4868 246.71520.40640.2342 0.5782-0
26、.4913-0.8521 -0.13240.33100.0340 0.6280R2=0.9815, F=0.442, p=0.0000 表2模型(2)的计算结果从表2得到最佳回归方程和用于统计分析的数据: (2)即模型98.15%可由变量解释,且|F|F(3,25)=0.004,p=0.000大于0.005,说明模型(2)整体是有效的,即国民收入、农业总产值、固定资产投资整体对财政影响显著5.2检查模型存在的序列相关性对于采用时间序列数据作样本,由于在不同样本点上解释变量以外的其他因素在时间上连续,带来它们对被解释变量的影响的连续性,所以往往存在序列相关性。模型(2)的拟合度已经达到了R2=0
27、.9815,但这个模型并没有考虑到我们的数据是一个时间序列。因此,为了排除随机项序列相关性对模型的影响,我们采用了残差分析法对模型(2)进一步改进。残差可以作为随机误差的估计值,画出etet-1的散点图(图9)能够从直观上判断的自相关性。残差数据见表3。年份1952195319541955195619571958残差-9.806719.05937.131519.3067-6.369619.7431-14.2544年份1960196119621963196419651966残差-2.7672-31.7878-3.178511.170821.4230 14.912433.6659年份1967196
28、819691970197119721973残差-27.2265-49.8054-3.49729.558722.959021.269521.2159年份1974197519761977197819791980残差-9.4577-27.1224-16.4919-1.911991.93428.3476-71.4846年份1952195319541955195619571958残差-12.73321.4793-5.70266.2455-21.60533.033411.4851年份1960196119621963196419651966残差-37.4752-1.1597-16.6442-7.478717
29、.223532.15615.043年份1967196819691970197119721973残差32.9908-22.6918-41.15255.414215.771428.712310.6054年份1974197519761977197819791980残差5.5024-24.5628-36.9943-14.23231.095186.33191.8876表3 模型(2)的残差图8 模型(2)etet-1的散点图从图8可以看出,大部分点落在第1,3象限,表明存在正得自相关。为了对的自相关性作定量诊断,并在确诊后得到新的结果,我们考虑如下的模型: , (3) 利用表3给出的残差,根据DW检验公
30、式计算出DW= 1.5106要根据DW的具体数值确定t是否存在自相关,则应该在给定的检验水平下,结合样本容量和回归变量数目,查D-W分布表,得到检验的临界值dL和dU,然后由DW所在的区间来决定。对于显著性水平=0.05,n=29,k=3,查D-W分布表,得到检验的临界值dL=1.27,和dU=1.56。正自相关0dL=1.27 dU=1.5624-dL 4-dU4无自相关负自相关不能确定不能确定DW=1.510666图10 与DW值对应的自相关状态现在DWdL,由图10中DW所在区间,我们可以看到DW落在不能确定的区间。我们分析原因可能是在1959-1962年,1966-1976年的重大历史
31、事件对几项因素,以及财政收入产生了特殊的影响,使数据的规律性降低,变得难以预测。但若是剔除这些年份,数据量将显得过少,以及会发生年代断代,因此我们依然保留原始数据。五、模型检验预测下面,我们将使用模型(2)对原始数据进行对比检验:年份1952195319541955195619571958实际值184216248254268286357预测值193.8067196.940240.868234.693274.369266.257371.254年份1959196019611962196319641965实际值444506271230266323393预测值477.109508.767302.788
32、233.178254.829301.577378.088年份1966196719681969197019711972实际值466352303447564638658预测值432.334379.226352.804450.497554.441615.041636.730年份1973197419751976197719781979实际值691655692657723922890预测值669.784664.458719.122673.491724.911830.066881.652年份1980实际值826预测值897.484 表4 财政收入的预测值与实际值对比可以看到,预测值与实际值还是相当吻合的。
33、故可以忽略模型中存在的时间序列相关性之后,我们查阅了1990年及2000年的国民收入、工业总产值、农业总产值、总人口、就业人口、固定资产投资,代入模型(2)。结果如下: 年份国民收入工业总产值农业总产值总人口就业人口固定资产投资199011412.9622452.27514611295463931.94517200098000.575710.6913873.612674372116.7732917.7年份财政收入(预测值)财政收入(实际值)19903539.682821.86200050961.713395.23 表5 1990年与2000年财政收入的实际值与预测值从表5可以看到,1990年预
34、测值与实际值相差大概33.35%,而2000年的预测值则偏差较大,高达67.12%。在查阅数据时,在固定资产投资一栏,我们查到了非常详细的分类,固定资产投资资金来源中国家预算内资金,固定资产投资资金来源中国内贷款,固定资产投资资金来源中自筹和其他资金等备注,并且随着社会经济的快速发展,人民币的币值也发生了极大的变化(这是可以来源于直观的生活经验的),模型对1990年和2000年预测值与实际值的匹配程度有限,但对于1952-1980年的经济数据来说,依然有一定的参考价值。六、模型的评价优点:1、 模型的决定系数较高,模型拟合效果非常好;2、 模型计算简单缺点:1、 采用面板数据建立线性回归模型存
35、在的多重共线性问题应用现在所有的技术都未能彻底解决,只能是剔除一些变量从而改善模型多重共线性,未能充分利用原有的数据;2、 模型自相关性未能进行有效的改善。3、 对相隔较远的时间预测效果吻合度低七、模型的推广与改进改进: 1、从前文的分析来看,我们认为1981年后,物价、通货膨胀、人民币币值等因素极大地影响了我们的预测数据,物价衡量指数,消费水平指数(CPI)等因素对国家财政收入的影响不可忽视,我们可以在模型中加入物价衡量指数,消费水平指数(CPI)因素的考虑,利用新的数据改善的拟合度。 2、我们要考虑全面不能轻易的删减任何数据,提高对数据的获取能力,使得模型的可信度更高。推广: 在我们对模型
36、的检验后发现,所建模型与财政收入实际值吻合度甚高,我们可以将模型应用到其他类似的经济问题中。八、参考文献1 刘卫国.MATLAB程序设计教程M.北京:中国水电水利出版社,2005.32 姜启源,谢金星,叶俊. 数学模型M.北京:高等教育出版社,2003.83 李子奈,潘文卿编著.-3版 计量经济学M.北京:高等教育出版社,2010.3 4CNKI中国统计年鉴数据库 九、附录 %部分程序1、描绘财政收入与各因素关系的拟合图MATLAB程序:fileimport datab2.xls %导入数据x1=b2(:,1); x2=b2(:,2); x3=b2(:,3); x4=b2(:,4); x5=b
37、2(:,5); x6=b2(:,6); y=b2(:,7); X=ones(30,1),x1,x2,x3,x4,x5,x6; X1=ones(30,1),x1; X2=ones(30,1),x2; X3=ones(30,1),x3; X4=ones(30,1),x4;X5=ones(30,1),x5;X6=ones(30,1),x6; 绘制拟合图:p1=polyfit(x1,y,1);xi1=500:3000;yi=polyval(p1,xi1);plot(x1,y,*,xi1,yi);title(国民收入与财政收入拟合图)figure;p2=polyfit(x2,y,1);xi2=0:700
38、0;yi=polyval(p2,xi2);plot(x2,y,*,xi2,yi);title(工业总产值与财政收入拟合图)figure;p3=polyfit(x3,y,1);xi3=400:1300;yi=polyval(p3,xi3);plot(x3,y,*,xi3,yi);title( 农业总产值与财政收入拟合图)figure;p4=polyfit(x4,y,1);xi4=55000:105000;yi=polyval(p4,xi4);plot(x4,y,*,xi4,yi);title(总人口与财政收入拟合图)figure;p5=polyfit(x5,y,1);xi5=20000:8000
39、0;yi=polyval(p5,xi5);plot(x5,y,*,xi5,yi);title( 就业人口与财政收入拟合图)figure;p6=polyfit(x6,y,1);xi6=0:600;yi=polyval(p6,xi6);plot(x6,y,*,xi6,yi);title(固定资产投资于财政收入拟合图)2、计算拟合相关数据:b1,bint,r,rint,stats=regress(y,X1) b2,bint,r,rint,stats=regress(y,X2) b3,bint,r,rint,stats=regress(y,X3) b4,bint,r,rint,stats=regres
40、s(y,X4) b5,bint,r,rint,stats=regress(y,X5) b6,bint,r,rint,stats=regress(y,X6)% b1返出是系数参值,stats返出第一个值是拟合系数的值,第二个值是MATLAB中置信度为0.05的F检验值,第三个值是相伴系数的值。b1 = 39.8438 0.3073bint = -8.0631 87.75070.2775 0.3371stats = 1.0e+003 *0.0009 0.4455 0 3.0848b2 = 202.1717 0.1106stats = 1.0e+003 *0.0009 0.3360 0 4.0129
41、b3 = -131.89770.8548stats = 1.0e+003 * 0.0008 0.1438 0.0000 8.5028b4 = -704.15770.0153stats = 1.0e+003 *0.0009 0.1971 0.0000 6.4877b5 = -70.96970.0173stats = 1.0e+004 *0.0001 0.0045 0.0000 2.0049b6 = 125.88851.3108stats = 1.0e+003 *0.0009 0.4643 0 2.9668disttool F= 0.004 %通过MATLAB软件中disttool工具箱查出F检验值F(1,29)=0.0043、fileimport datab1.xlsx1=b1(:,1);x2=b1(:,2); x3=b1(:,3); x4=b1(:,4);