kevinj0419

Untitled

Jan 5th, 2016
83
0
Never
Not a member of Pastebin yet? Sign Up, it unlocks many cool features!
MatLab 6.90 KB | None | 0 0
  1. clear;clc;
  2.  
  3. file='20160105_B-N_NR_X';
  4.  
  5. q=0.9;
  6. a=0.68;
  7. D=0.075;
  8. c0=0.1;
  9. c1=0.001;
  10. X=2^(0.5);
  11.  
  12. for i=1:1:101
  13.     for j=1:1:101
  14.        
  15.         b=0.01*i-0.01;
  16.         G=(j-j)/100;
  17.        
  18.         pB1 = (q*(a + b))/2 - (q*(q*(a + b) - 2*c0))/(8*(q - c1));
  19.         wB1 = (q*(2*c0 + a*q + b*q - 2*c1*(a + b)))/(4*(q - c1));
  20.         cB1 = ((4*c0 - a*c1 - b*c1)*q - 2*c0*c1)/(4*(q - c1));
  21.             MSB_a1 = (q*(a + b) - 2*c0 + 4*(a - b)*(q - c1))/(8*(q - c1));
  22.             MSB_b1 = (q*(a + b) - 2*c0 - 4*(a - b)*(q - c1))/(8*(q - c1));
  23.             MSB1 = ((a + b)*q - 2*c0)/(4*(q - c1));
  24.        
  25.         wN1 = -(q*(q*(a + b) - 2*c0 - 4*(q - c1)*(a + b)))/(8*(q - c1)) - (q*(q*(a + b) - 2*c0)^2)/(4*(q - c1)*(q*(a + b) - 2*c0 + 4*(q - c1)*(a - b + 2*D)));
  26.         pmN1 = (q*(a + b))/2 + (q*(2*c0 - 4*D*q - 3*a*q + b*q))/(8*(q - c1)) - (q^2*(q*(a + b) - 2*c0))/(8*(q - c1)^2);
  27.         cN1 = (4*c0*c1^2 + 8*c0*q^2 - 10*c0*c1*q - 4*D*c1*q^2 + 4*D*c1^2*q - 5*a*c1*q^2 + 4*a*c1^2*q - b*c1*q^2)/(8*(q - c1)^2);
  28.             MSN_a1 = (q*(a + b) - 2*c0 + 4*(a - b)*(q - c1))/(8*(q - c1));
  29.             MSN_b1 = b/2 - a/2 + (2*c0 - 4*D*q - 3*a*q + b*q)/(8*c1 - 8*q) - (q*(2*c0 - q*(a + b)))/(8*(c1 - q)^2);
  30.             MSN1 = (q*(a*q - 2*c0 + b*q))/(8*(c1 - q)^2) - ((D*q)/2 - c0/2 + (a*q)/2)/(c1 - q);
  31.        
  32.         pB2 = (q*(a + b))/2 - (q*(q*(a + b) - 2*c0))/(8*(q - c1));
  33.         wB2 = (q*(2*c0 + a*q + b*q - 2*c1*(a + b)))/(4*(q - c1));
  34.         cB2 = ((4*c0 - a*c1 - b*c1)*q - 2*c0*c1)/(4*(q - c1));
  35.             MSB_a2 = (q*(a + b) - 2*c0 + 4*(a - b)*(q - c1))/(8*(q - c1));
  36.             MSB_b2 = (q*(a + b) - 2*c0 - 4*(a - b)*(q - c1))/(8*(q - c1));
  37.             MSB2 = ((a + b)*q - 2*c0)/(4*(q - c1));
  38.    
  39.         wN2 = (q*(a + b))/2 + (q*(a*q - 2*c0 + b*q))/(4*(c1 - q));
  40.         pmN2 = (q*(a + b))/2 + (q*(a*q - 2*c0 + b*q))/(8*(c1 - q)) - (3*q^2*(a*q - 2*c0 + b*q))/(16*(c1 - q)^2);
  41.         cN2 = (8*c0*c1^2 + 16*c0*q^2 - 18*c0*c1*q - 7*a*c1*q^2 + 4*a*c1^2*q - 7*b*c1*q^2 + 4*b*c1^2*q)/(16*(c1 - q)^2);
  42.             MSN_a2 = a/2 - b/2 - ((a*q)/8 - c0/4 + (b*q)/8)/(c1 - q);
  43.             MSN_b2 = b/2 - a/2 - ((a*q)/8 - c0/4 + (b*q)/8)/(c1 - q) + (3*q*(a*q - 2*c0 + b*q))/(16*(c1 - q)^2);
  44.             MSN2 = -((4*c1 - 7*q)*(a*q - 2*c0 + b*q))/(16*(c1 - q)^2);
  45.  
  46.            
  47.         if b >= (D + pB1/q) & pB1 >= 0 & wB1 >= 0 & cB1 >= 0 & MSB_a1 >= 0 & MSB_b1 >= 0 & pB1 - wB1 >= 0 & wB1 - cB1 >= 0
  48.             B1 = (q*(a + b) - 2*c0)^2/(16*(q - c1));
  49.         else
  50.             B1 = 0;
  51.         end
  52.        
  53.         if b >= (D + pB1/q) & pB1 >= 0 & wN1 >= 0 & pmN1 >= 0 & cN1 >= 0 & MSN_a1 >= 0 & MSN_b1 >= 0 & pB1 - wN1 >= 0 & wN1 - cN1 >= 0 & pmN1 - cN1 >= 0
  54.             N1 = ((16*D^2 + 24*D*a - 8*D*b + 13*a^2 + 2*a*b + 5*b^2)*q^4 + (- 32*c1*D^2 - 40*c1*D*a + 24*c1*D*b - 16*c0*D - 20*c1*a^2 - 28*c0*a - 12*c1*b^2 - 12*c0*b)*q^3 + (16*D^2*c1^2 + 16*D*a*c1^2 - 16*D*b*c1^2 + 16*D*c0*c1 + 8*a^2*c1^2 + 40*a*c0*c1 + 8*b^2*c1^2 + 24*b*c0*c1 + 20*c0^2)*q^2 + (- 32*c0^2*c1 - 16*a*c0*c1^2 - 16*b*c0*c1^2)*q + 16*c0^2*c1^2)/(- 64*c1^3 + 192*c1^2*q - 192*c1*q^2 + 64*q^3);
  55.         else
  56.             N1 = 0;
  57.         end
  58.        
  59.         if b < (D + pB1/q) & b >= pB1/q & pB1 >= 0 & wB1 >= 0 & cB1 >= 0 & MSB_a2 >= 0 & MSB_b2 >= 0 & pB1 - wB1 >= 0 & wB1 - cB1 >= 0
  60.             B2 = (q*(a + b) - 2*c0)^2/(16*(q - c1));
  61.         else
  62.             B2 = 0;
  63.         end
  64.  
  65.         if b < (D + pB1/q) & b >= pB1/q & pB1 >= 0 & wN2 >= 0 & pmN2 >= 0 & cN2 >= 0 & MSN_a2 >= 0 & MSN_b2 >= 0 & pB1 - wN2 >= 0 & wN2 - cN2 >= 0 & pmN2 - cN2 >= 0
  66.             N2 = -((a*q - 2*c0 + b*q)^2*(16*c1^2 - 32*c1*q + 25*q^2))/(256*(c1 - q)^3);
  67.         else
  68.             N2 = 0;
  69.         end
  70.  
  71.        
  72.         Y = [B1 N1 B2 N2];
  73.         Z = max(Y);
  74.         if Z == 0
  75.            
  76.             AA=[{'N/A'},NaN,NaN,NaN,NaN,NaN,NaN,NaN,NaN,NaN,NaN,NaN,NaN,NaN,NaN,NaN,NaN,NaN,NaN,NaN,NaN,NaN,NaN,NaN,NaN,NaN,NaN,NaN,NaN,NaN,NaN,NaN,NaN,NaN,NaN,NaN,NaN,NaN,NaN,NaN,NaN,NaN,NaN,NaN,NaN];
  77.             str = {'N/A'};
  78.             hold on
  79.             plot(b,a,'.w');
  80.         else
  81.             if Z == B1
  82.                 AA=[{'B1'},pB1,wB1,cB1,pB1-wB1,wB1-cB1,wN1,pmN1,cN1,pB1-wN1,wN1-cN1,pmN1-cN1,NaN,NaN,NaN,NaN,NaN,NaN,NaN,pB2,wB2,cB2,pB2-wB2,wB2-cB2,wN2,pmN2,cN2,pB2-wN2,wN2-cN2,pmN2-cN2,MSB_a1,MSB_b1,MSB1,MSN_a1,MSN_b1,MSN1,NaN,NaN,NaN,MSB_a2,MSB_b2,MSB2,MSN_a2,MSN_b2,MSN2];
  83.                 str = {'B1'};
  84.                 hold on
  85.                 plot(b,MSB_a1,'.r');
  86.                 plot(b,MSB_b1,'.g');
  87.                 plot(b,MSB1,'.b');
  88.             end
  89.             if Z == N1
  90.    
  91.                 AA=[{'N1'},pB1,wB1,cB1,pB1-wB1,wB1-cB1,wN1,pmN1,cN1,pB1-wN1,wN1-cN1,pmN1-cN1,NaN,NaN,NaN,NaN,NaN,NaN,NaN,pB2,wB2,cB2,pB2-wB2,wB2-cB2,wN2,pmN2,cN2,pB2-wN2,wN2-cN2,pmN2-cN2,MSB_a1,MSB_b1,MSB1,MSN_a1,MSN_b1,MSN1,NaN,NaN,NaN,MSB_a2,MSB_b2,MSB2,MSN_a2,MSN_b2,MSN2];
  92.                 str = {'N1'};
  93.                 hold on
  94.                 plot(b,MSN_a1,'.r');
  95.                 plot(b,MSN_b1,'.g');
  96.                 plot(b,MSN1,'.b');
  97.             end
  98.             if Z == B2
  99.  
  100.                 AA=[{'B2'},pB1,wB1,cB1,pB1-wB1,wB1-cB1,wN1,pmN1,cN1,pB1-wN1,wN1-cN1,pmN1-cN1,NaN,NaN,NaN,NaN,NaN,NaN,NaN,pB2,wB2,cB2,pB2-wB2,wB2-cB2,wN2,pmN2,cN2,pB2-wN2,wN2-cN2,pmN2-cN2,MSB_a1,MSB_b1,MSB1,MSN_a1,MSN_b1,MSN1,NaN,NaN,NaN,MSB_a2,MSB_b2,MSB2,MSN_a2,MSN_b2,MSN2];
  101.                 str = {'B2'};
  102.                 hold on
  103.                 plot(b,MSB_a2,'.r');
  104.                 plot(b,MSB_b2,'.g');
  105.                 plot(b,MSB2,'.b');
  106.             end    
  107.             if Z == N2
  108.  
  109.                 AA=[{'N2'},pB1,wB1,cB1,pB1-wB1,wB1-cB1,wN1,pmN1,cN1,pB1-wN1,wN1-cN1,pmN1-cN1,NaN,NaN,NaN,NaN,NaN,NaN,NaN,pB2,wB2,cB2,pB2-wB2,wB2-cB2,wN2,pmN2,cN2,pB2-wN2,wN2-cN2,pmN2-cN2,MSB_a1,MSB_b1,MSB1,MSN_a1,MSN_b1,MSN1,NaN,NaN,NaN,MSB_a2,MSB_b2,MSB2,MSN_a2,MSN_b2,MSN2];
  110.                 str = {'N2'};
  111.                 hold on
  112.                 plot(b,MSN_a2,'.r');
  113.                 plot(b,MSN_b2,'.g');
  114.                 plot(b,MSN2,'.b');    
  115.  
  116.                 xlabel('Beta');
  117.                 ylabel('MS')
  118.             end
  119.         end
  120.         matrix((i-1)*101+j,:)=[q,a,c0,c1,b,D,str,AA(2),AA(3),AA(4),AA(5),AA(6),AA(7),AA(8),AA(9),AA(10),AA(11),AA(12),AA(13),AA(14),AA(15),AA(16),AA(17),AA(18),AA(19),AA(20)...
  121.             ,AA(21),AA(22),AA(23),AA(24),AA(25),AA(26),AA(27),AA(28),AA(29),AA(30),AA(31),AA(32),AA(33),AA(34),AA(35),AA(36),AA(37),AA(38),AA(39),AA(40),AA(41),AA(42),AA(43),AA(44),AA(45)];
  122.     end
  123.    
  124. end
  125.  title={'q','Alpha','c0','c1','Beta','Delta','OEM Strategy','pB1','wB1','cB1','pB1-wB1','wB1-cB1','wN1','pmN1','cN1','pB1-wN1','wN1-cN1','pmN1-cN1','pF1','wF1','pmF1','cF1','pF1-wF1','wF1-cF1','pmF1-cF1'...
  126.      ,'pB2','wB2','cB2','pB2-wB2','wB2-cB2','wN2','pmN2','cN2','pB2-wN2','wN2-cN2','pmN2-cN2','MSB_a1','MSB_b1','MSB1','MSN_a1','MSN_b1','MSN1','MSF_a1','MSF_b1','MSF1','MSB_a2','MSB_b2','MSB2','MSN_a2','MSN_b2','MSN2'};
  127.  sheet=1;
  128.  xlRange = 'A2';
  129.  xlswrite(file,title,sheet,'A1');
  130.  xlswrite(file,matrix,sheet,xlRange);
Add Comment
Please, Sign In to add comment